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ABSTRACT 



In this thesis, an adaptive lattice algorithm is derived for an AR\1A digital lattice 
filter, whose parameters are estimated using a generalized .Ylullis- Roberts criterion for 
parameter estimation. Design of the ARMA lattice filter based on this generalized cri- 
terion is studied as is the accuracy of the parameter estimation algorithm used in its de- 
sign. Application of the derived lattice algorithm to system identification modeling is 
demonstrated through computer simulation of various system identification problems. 
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I. INTRODUCTION 



Digital signal processing is a field which is rapidly expanding due to advances in 
modern technology. Essential to this field are digital filters. Modeling these filters con- 
stitutes much of the effort involved in digital signal processing. The filters provide a 
transfer function which describes the relationship between filter input and output. 
Autoregressive (AR), moving average (MA) and combination autoregressive moving 
average (ARMA) models are widely used to represent the transfer function of a digital 
filter. A filter transfer function is commonly described in direct form. This form is a ratio 
of two polynomials, usually of the form, 



//(-) 



Biz) 

B(z) 



ftp +_Bj_z ' 4- hj\z 
1 4- z + ••• + a s z 



(I- 1) 



The above equation describes an ARMA model of order (s,t) where s is the order of the 
denominator and t the order of the numerator. The a s parameters form the 
autoregressive portion of the ARMA model. The b, parameters form the moving average 
portion of the ARMA model. If all the autoregressive parameters are zero, then the filter 
transfer function H(z) is strictly a moving average process of order t. If all the moving 
average parameters are zero except for b 0 equal to one, then the filter transfer function 
is strictly an autoregressive process of order s. 

A. OBJECTIVES OF THE THESIS 

Fundamental to the design of digital filters is estimation of AR, MA or ARMA 
parameters. Accurate and efficient parameter estimation has been the subject in much 
of the related digital signal processing literature [Refs. 1 ,2.3]. The first objective of this 
thesis is to confirm the proposed ARMA parameter estimation algorithm of [Ref. 4: pp. 
619-621], which leads to the design of a new ARMA digital lattice filter. The proposed 
algorithm is a generalization of the Mullis-Roberts criterion for parameter estimation 
known as the modified least squares problem [Ref. 5: pp. 227-22S]. The algorithm uses 
two recursive formula to estimate the parameters. One is an A R recursive formula which 
estimates ARMA parametes as the AR order is increased by one. The other is an MA 
recursive formula which estimates ARMA parameters as the MA order is increased by 
one. This algorithm is unique in that it allows for the design of an ARMA model with 



1 



arbitrary AR and MA orders with no dependency of an AR model on an MA model or 
vice versa. The AR.MA digital filter designed from the proposed ARMA parameter es- 
timation algorithm is in the form of a lattice structure. Lattice realizations of filters are 
widely used and provide excellent analysis of prediction errors [Refs. 6 : pp. 165-168.7). 
Gray and Markel developed an algorithm which produces lattice realizations of filters 
from the direct form [Ref. 8). 

The second objective of the thesis is to make the proposed ARMA digital lattice 
filter of [Ref. 4: p. 662J adaptive. An adaptive lattice filter is one in which the lattice 
coefficients are automatically adjusted by an adaptive algorithm to yield the optimum 
filter design. The adaptive lattice algorithm derived in this thesis is based on the widely 
used least mean square (LMS) algorithm. Adaptive filters have many applications [Ref. 
9: pp. 7-31] including. 

1. System identification. 

2. Digital representation of speech. 

3. Adaptive auotoregressive spectrum analysis. 

4. Adaptive detection of a signal in noise of unknown statistics. 

5. Echo cancellation. 

6. Adaptive line enhancment. 

7. Adaptive beamforming. 

The need for an adaptive filter is made apparent by considering a filter of fixed design 
which is optimized for given input conditions. In practice, the complete range of input 
conditions may not be known or could change from time to time. A filter of fixed design 
would not produce optimum results under these conditions. An adaptive filter, which 
yields optimum results given changing input conditions, will give superior performance 
to one of fixed design. 

The last objective is to analyze convergence properties of the derived adaptive lattice 
algorithm. This is accomplished by computer implementation of the adaptive algorithm. 
The output of a known transfer function is compared to the output of the adaptive 
lattice filter given a common input. Plots of the error between the two outputs and 
lattice coefficient convergence are obtained. 

B. ORGANIZATION OF THE THESIS 

Chapter II is designed to present the ARMA parameter estimation algorithm and 
ARMA digital lattice filter proposed in [Ref. 4: pp. 617-628). Computer simulation of the 



algorithm was performed and results are shown. A brief review of the Mullis- Roberts 
criterion is provided to establish a reference for expanding this criterion to the proposed 
ARMA parameter estimation algorithm. An adaptive lattice algorithm is derived in 
Chapter 111 which makes the proposed ARMA digital lattice filter adaptive. The adap- 
tive lattice algorithm is efficient and accurate. Chapter IV contains experimental results 
which show convergence aspects of the adaptive lattice algorithm when applied to 
ARMA lattice filters. Conclusions about the proposed ARMA parameter estimation 
algorithm as well as the derived adaptive algorithm are discussed. 
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II. ARMA DIGITAL LATTICE FILTER 



In this chapter, we will review the Mullis-Roberts criterion for solving linear ap- 
proximation problems and introduce analysis equations of the ARMA digital lattice fil- 
ter. The criterion used in the formulation of the ARMA digital lattice filter is a 
generalized form of the Mullis-Roberts criterion [Ref. 5: pp. 227-228], which has been 
given as a modified least mean square problem for ARMA parameter estimation. 

A. MULLIS-ROBERTS CRITERION 

The Mullis-Roberts criterion evolved from considering second order statistics in 
conjuction with first order information about a process to obtain filter approximations. 
Consider the bounded impulse response sequence h = {/? 0 , h t , ... } containing first order 
information about the filter h having a frequency response function, 



h k e~ Jcok (2.1) 

k = 0 

Let ( u, } be a zero-mean, unit-variance, white-noise sequence and ( y, } be the output 
process corresponding to the input } u, } . then we have the following convolutional re- 
lationship given by, 



oo 

y I = Y J h k“ l -k ( 2 - 2 ) 

A-=0 

Second order information about the filter h is obtained from the autocorrelation se- 
quence { r k } given as 



>'k 



E tv t y,+k) = YjW 



■k+i 



1=0 



(2.3) 



From equations (2.1) and (2.2), the second order interpolation problem is to find a 
lowest order recursive filter which matches the data { /i 0 . ... , h m , r 0 r m } , where h, re- 

presents the first order information and r, the second order statistics. 
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Let us now consider the case where only first order information about a process is 
known. That is. given an impulse response sequence { h 0 , ft ,. ... } , we want to find a re- 
cursive filter of the form, 




% + 7U 1-1 f %£ 

1 + UjZ + ••• + a n z 



(2.4) 



which approximates H(z) and therefore the impulse response sequence { /i 0 . h u ... } . We 
also desire the frequency response H{e ;al ) to approximate the desired response II(e JW ). 
Suppose that H(e jw ) is chosen such that it minimizes the integral squared error, 



2 7c J 



\Ilie J( °) - H(e‘ 



Ja> 



fdco = ||/1-/1|| 2 



(2.5) 



Using the Parseval relation, we can obtain an alternative definition of the approximation 
error. 



2- J 



1/7(0 - lli^tfdu = |!/i-/i|| 2 = ) (h k -h k f 



(2.6) 



k=0 

If the filters Il{z) and H{z) are driven by the same white noise source, equation (2.6) 
describes the output error between the filters which we write as, 

P-Zill 2 = EiS'-vf = E(ef (2.7) 



where y, and y, are the outputs of the respective filters when driven by the same white 
noise source as in (2.2). Minimizing (2.7) is a nonlinear programming problem requiring 
the entire impulse response sequence. As a result, computational efforts for obtaining a 
solution are inefficient. A modification to the problem was introduced [Ref. 5: pp. 
227-228] which considered a cost function that is quadratic in the coefficients of the re- 
cursive filter given by (2.4). The modification is described as follows. Let 

<?(*) = ‘ A (7o + (2-8) 
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and 



A{ 2 ) = r V ( 1 + a,.-" 1 + - + a n z~ n ) (2.9) 

be the numerator and denominator polynomials, respectively, in (2.4), where ,V= max 
(m.i ;). The task now is to find coeflicients which minimize the quadratic form, 

E(c ,) 2 = 3T J“ \II{t J<u )A(e fu> ) - (2. 10) 

This is a standard quadratic minimization problem whose integral can be expressed in 
terms of the coefficients of polynomials A(z) and Q(z) and the data 

{/i# h m , r 0 . ... , r m ). relating to the filter 7/(r). This problem is shown in Figure 1 and 

equation (2.10) is known as the Mullis-Roberts criterion. 



u t 


H(z) 




A(z) 










▼ e 












+ t 

















Q(z) 





Figure I. Modified least squares problem 



B. GENERALIZED MULLIS-ROBERTS CRITERION 

In order to define the new criterion used for ARM A parameter estimation we con- 
sider the following transfer function with input sequence ( x(k), k = 1,2, ... } and output 
sequence { y(k), k = 1,2, ... } written as, 



m 



H v (z) 

//,(-) 



( 2 . 11 ) 
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Figure 2. Equivalent input/output model 



where H>{ r) and H,(z) are reference polynomials which we desire to model. An equivalent 
model is shown in Figure 2 where u(k) is an intermediate signal, and the realization is 
similar to that of the direct form realization 11 [Ref. 10: p. 151]. Let the intermediate 
sequence { u(k). k = 1,2, ... } be a zero-mean white gaussian process. The model of Fig- 
ure 2 can then be transformed into the model of Figure 3 with u(k) as the common input 
to both transfer functions II y ( :) and //.(r). 



u (k) 






y(k) 




Hy (2) 






t 








x(k) 






H x (z) 

















Figure 3. Transformed model 



Note that we earlier defined transfer functions //,(r) and H t {z) as polynomials of the 
reference model. For each transfer function an estimation polynomial is defined such 
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that 



A{:) — 1 4- t7]C ' + ••• + a s z 


estimates Il x (z) 


(2.12) 


Biz) = b 0 + V " 1 + - + b,z" 


estimates // v (c) 


(2.12) 



where a, and b„ (/= 1.2, ... , s) and (j = 1,2 /) , are the AR and MA parameters, re- 

spectively. of the combined ARM A model formed by A(z) and B(.z). This refined least 
squares problem is shown in Figure 4 and is a generalized form of the Mullis-Roberts 
criterion. If the reference model polynomial HJy) in Figure 4 is equal to unity, we obtain 
the Mullis-Roberts criterion shown in Figure 1. Therefore, by including reference 
polynomial //,(;). the new criterion for ARMA parameter estimation becomes, 



= //,(:)5(:)| : T (2.14) 
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Figure 4. Refined least squares problem 



Minimizing £ M is accomplished by calculating the coefficients of A(z) and B(z) which 
minimize e(k) : of Tigure 4. Another form of eq (2.14) is obtained by applying Parseval's 
theorem and is expressed as, 

E SJ = E[(A(z)y(k) - B(z)x(k)) 2 ] (2.15) 

which is obvious from Figure 4. The coefficients a„ (/ = 1,... . 5 ), and b r (j = 1 r), 

which minimize can then be calculated using the normal equations for ARMA 
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parameter estimation. In order to obtain the normal equations for the problem in (2.15). 
let us define the following: 



a w = [>i ••• oj and b s<t = [6, ... Z> f ] (2.16) 

are the vectors of AR and MA parameters, respectively, 

M*) = [ y(k) ...y(k - s) -x(k ) ... -x(k - t) ] (2.17) 

is the data vector consisting of both input and output data elements and 

R s ,= Elh s>t (k) r h s<l (k)-} (2.18) 

is the data autocorrelation matrix. The criterion is to minimize the mean squared error 
E s , t = E [ e\k) ] = £[[[ 1 a s t b 0 b,, r ] h£ ] : ] 

L (~> 19) 

= e[ [>(*) + [a w b Q b SJ ] hj, ] 2 ] 

Now following the standard calculus of variation optimization procedure for minimizing 
£)„ [Ref. 11: pp. 100-110], yields the normal equations in the matrix form, 

[ 1 a J>f b 0 b w ]R v = [ min E sr 0 0 0 ] (2.20) 

It is interesting to note that if E SJ in equation (2.14) is zero, then we have the following 
equality, 



H(z) 



IE(z) 

//,(-) 



B(z) 

A(z) 



( 2 . 21 ) 



so the estimate for the total reference model H(z) is the ratio of B(z), the \1A part and 
A(z), the AR part of an estimated ARMA model. 

C. ARMA PARAMETER ESTIMATION 

Now that the criterion for ARMA parameter estimation has been established, the 
solution method to estimate the ARMA model parameters to minimize equation (2.14) 
is considered. Let x(k) and y(k) be the input and output signals, respectively, of the es- 
timated ARMA model. Using a difference equation representation, this process is de- 
scribed by 
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( 2 . 22 ) 



S t 

y(k) = -Vfl ; y(k -j) + Y b, x(k - i) 
j= 1 /=o 

For these input and output signals we define four estimation, or prediction, models as 
follows. The forward estimation signal for ;c(A) is defined as, 

t s 

Xjik) = ~Y b ? x(k -0 + Y a J y {k -j) (2.23) 

;=i 7=i 

where b’ and a* are the corresponding estimation parameters. The forward esitmation 
signal forv(A) is similarly defined as. 






)(k) = - V aJ y(k -j) + V b y x{k _ /) 

The backward estimation errors for x(k — /) and y(k — s) are then given by, 

f-I s - 1 

4(* - /) = - y Z>f x(k - i) + y of y(A -J) 



(2.24) 



(=0 



1-1 



y=o 



/-i 



j*(^ - *) = - 2^ a j y( k - j) + Y x ^ k ~ 

j = o 



(2.25) 



(2.26) 



/=o 



where the superscripts g and d indicate the backward estimation parameters for .v and 
y, respectively. 

From these estimation models, we can now obtain the four prediction errors at any 
given time k. These errors are expressed as differences between the predicted value and 
actual value of an input or output signal, namely, 



&*) - -*(*) + ty*) 


(2.27) 


y? 

ii 

V. 

i 


(2.2S) 


ii 

y? 

i 

i 

i 


(2.29) 
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) 



d s ,'(k)=y(k-s)-y b (k-s) 



(2.30) 



We now use the vector notation to simplify the expressions for prediction errors. In the 
following, the forward error elements corresponding to both .v and y form a vector 
x S ',{k), given by 

V(*) = [ -</*) v-;;,(/0] = h v (A-) c£ (2.31) 

and the backward error vectors are given by 

Ss/k) = - hjk) Gj, (2.32) 

d s , t (k) = hjk) D l, (2.33) 

where C u , and G Jif and D ;[ are the forward estimation parameter matrix and backward 
estimation parameter vectors, respectively, defined as 







~0 


X 

<7, ... 


a , 


1 


b* 


... bf 
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V 

a \ - 
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... *J 
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.34) 


G w = 


Oo° 


a s 

... a i _ ] 


i 0 


f ° 

bo 




» ° 

u t - 1 


a 


(2. 


.35) 


= 


1 — 1 


d 

... 


1 


i d 

bo 




id 
"t - 1 


0] 


(2. 


.36) 



It can be shown that the prediction errors satisfy some orthogonal conditions. These 
conditions are similar to those found in AR modeling problems [Ref. 6 : pp. 116-121]. 
We now list the orthogonality conditions for the ARMA formulation in discussion 
without proof as the following: 

e[ »£(*) y(k-j) ] =0 £[ vi'(k) y(k-j) ] =0 

E[ v*(£) *(*-/) ]=0 £[ ii(k) x(k-i) ] =0 

££&,(*-!) }<k -j)J = 0 (3.37) 



£[&,(*-!) 4k ~ /)] = 0 



Eld u {k- 1) y(k -j)-] = 0 
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Eld s /k- 1) x(k - /)] = 0 



where i = 1,2. ... , t and j = 1,2, ... , 5. 

In Section B we have obtained a set of normal equations in terms of R,, and the 
AR.MA estimation parameters. In this section, we have defined four sets of forward and 
backward estimation parameters and established some orthogonal relationships. In what 
follows, we derive a set of equations which relates the coefficients of the estimated 
AR.MA model with those of the forward estimation parameter matrix C I( . Consider the 
expected value of the forward prediction error, \„(k) and the data h Si ,(A). Since the pre- 
diction error is orthogonal to all past samples of data y(k —j), x(k — i ) but not to y(k) 
or x(k) as listed in (2.37), the result is a matrix which is defined as 



E[y s , t (k) T h sl (k) ] = 




0 c 2 0 
0 c. 0 



(2.3S) 



where 

c, = -E [ v?/k)y(k) ] = -E[ £(/<) <,(/c) ] = E*J (2.39) 

is the crosscorrelation between the forward prediction errors of .r and y at lag zero. 

c 2 = £ [ v*,(k) x(k) ] = £ [ {'•s/k)) 2 ] = Efj (2.40) 

is the forward prediction error power of .r(A) 

c 3 = E[xi'(k)y(k) ] = £[ « f (A)) 2 ] = El, (2.41) 

is the forward prediction error power of y(k) and 

L = -E[ V £,(*) X(k) ] = -£ [ vl,(k) vfjik) ] = E?f (2.42) 

is the crosscorrelation between the forward prediction errors of v and .v at lag zero. In 
another interpretation, the left hand side of equation (2.3S) can be written as, 

E [ VsJLk) T h s Jik) ] = £ [ C w hjkfhjk) ] = C w R v (2.43) 



and we have 
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(-5,! ^5,: — 



(2.44) 



£?/ o e*, o 

E h 0 E U o 



from equation (2.38). 

A similar approach for both backward prediction errors yields the following 
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T7S d 
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F° 
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pgd 

Aj.r 



(2.45) 



In order to express the coefficients of the ARMA estimation model in terms of the co- 
efficients of the parameter estimation matrix C s> , and parameter vectors G Si , and D s , , 
consider the combination of the normal equation (2.20) and the parameter estimation 
matrix and vectors. From equation (2.44) the normal equation for ARMA parameter 
estimation may be written as 



0 a v 1 b£, 

1 a ),r ^ ? 



R 



S,l 



E U 0 E* t 0 

E h 0 E*J 0 



(2.46) 



Combining equation (2.46) with (2.20) we obtain 



" 1 
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*0 


b,/ 




£ if min 
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E xy 
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E x 
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0 


E xy 

^ s,t 
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E*y 

In equation (2.47) multiply row 2 by — — then subtract row 2 from row 3 and equate 
the result to row 1. We then obtain the following relations 



f xy 

J 'S.t 



_X 

a 5,f r-JC 

E s,t 


(2.4S) 


E xy 

E x 

~ s,t 


(2.49) 


E xy 

b i,/ V X 
Ee t 


(2.50) 



'V 



and, 
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(2.51) 



E s , min = E y s l 




Calculation of the ARM A estimation model coefficients requires knowledge of the esti- 



calculate the estimation parameters as the AR or MA order of the estimation model in- 
creases by one. Given the parameters of the ARM A estimation model of order ( 5 , /) 
these formulas calculate the (s + 1 , /) or (s , l + 1) ARMA parameters. For a compre- 
hensive derivation of these recursive formulae see [Ref. 4: pp. 619-621]. The AR-type 
recursive formula for the forward parameter estimation matrix and backward estimation 
parameter vectors is 



The matrices I, and I 2 are of dimension (s + r + 2 x s + / + 3). They are introduced to 
provide symmetry to the matrix algebra and preserve initial condition calculations. We 
design the matrices to perform the following operations 



mation parameters a’ r , aj >f , bf„ b? r , b 0 and the values £* r , £>,. £’> . Recursive formulas 




(2.52) 



where 



u, = - (4r‘ ct, t 2 ] 




(2.53) 



and the (s+ l,t) prediction error powers are recursively calulated as follows 




(2.54) 
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C w, 


... (L ) J+1 C0 J+ 2 ••• ^ ^ 


1 = [ CO, .. 


• <^+i ^ W 5+2 


(2.55) 


[ CO, . 


.. CJ 5+1 C£) j _ ) _2 ••• <^+7+2 3 1 2 


II 

1 — 1 
o 

e 


... 0 CO J+ 2 ^Vl-r+l 1 


(2.56) 



The values t, through t 4 can be obtained through the correlation data and forward and 
backward estimation parameters. We express them mathematically as 

C tj t 2 3 = £ [ d s l {k — 1) v i( ,(A) ] 

t 3 = -£[>(* -s-1)^, (A) ] (2.57) 

t 4 = £ C y{k - s - 1) d s ,{k) ] 

The MA-type recursive formula for the forward estimation parameter matrix and back- 
ward estimation parameter vectors is obtained in a similar manner. The recursive for- 
mula is given by 

C j,/+1 = G s,t *3 + 11 T G j.f 

d 4 ,/+i = [ D i,/ + G s,t ] h (2.58) 

G s,t+ 1 = G S .1 + [ n 3 G s,i + , + n 5 G s , ] I 3 



where 



■*l - - (£■*)-’ [ T-, r' ; ] 



sd 



)U = 



-I - 



Eh 

n 3 = - [ T 'l T> 2 ] E v 

(Efj v' 4 — E s s l t' 3 ) 



n A = • 



[ (£f /) 2 - 4 ., 4 ] 

n 5 = n A n 2 



(2.59) 



and the (s,t+ 1) prediction error powers are calculated using the following recursive for- 
mulas 



E i,/+i = E v + n rO'i z ' 2 ] 

E ?/ +] = t ' 3 + « 2 t ' 4 

Ei ,+ 1 = El, + [ t', t' 2 ] n 3 r + » 4 t'3 + « 5 t' 4 
4+> = El, + w 2 Elf 



(2.60) 
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where I 3 and I 4 are (s + / -f 2 x s + / + 3) dimensional matrices which we have designed 
to perform the following operations 

[ w, ... <u J+1 co, +2 ... co s+l+2 2 I 3 = [ - Wj+i - ^s+t+2 0 ] (2-61) 

[ co, ... co J+1 co, +2 ... co i+f+2 J I 4 = [ 0 co, ... co, 0 c o, +2 ... co i+f+2 ] (2.62) 

The values of t', through x\ are calculated using correlation data in conjunction with 
current forward and backward parameter estimation values. We express these quantities 
mathematically as 

[t'i T' 2 ] = -E[g v (A-l)v v (A)] 

t' 3 = -£[.v(A-/ - 1) d SJ (k) ] (2.63) 

r'i= E [ x(k - t - 1) g Sil {k) ] 

It is interesting to note that the \1 A- type recursive formula is the complimentary form 
of the AR-type formula and that the two are identical if the variables associated with the 
input signal ,v(A) and the variables associated with the ouput signal y(k) are interchanged. 
That is. we replace y{k). G J( and g s Jk) "’ith —x (k), D 5i , and d St ,{k) and vice versa. 

1. Experimental Results 

The ARMA parameter estimation algorithm of [Ref. 4: pp. 619-621] based on 
the recursive formulas of equations (2.52) and (2.5S) was implemented using the Fortran 
program found in Appendix A. This program calls subroutines which compute the 
ARMA model parameters as the AR order is increased by one and as the MA order is 
increased by one. These subroutines are shown in Appendix B and Appendix C, respec- 
tively. In the main program, an input data sequence of white Gaussian noise is passed 
through a known reference model producing an ouput data sequence. We obtain 
autocorrelation and crosscorrelation data from these input and output sequences. The 
correlation data is used to calculate initial values of the error powers for .v and y as well 
as r, through r 4 and t', through r' 4 . Next we obtain estimates of the reference model by 
employing the recursive formulas (2.4S) through (2.50), (2.52) and (2.5S). Several refer- 
ence models were estimated beginning with a strictly AR process of order s = 4 having 
as its transfer function 



1 — 0.2 z _I + 0.62 z~ 2 — 0.152 z~ 3 + 0.3016 z - " 1 
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The actual values of the reference model parameters and the ARMA model parameters 
which estimate this reference model are listed below 





ACTUAL 


ESTIMATED 


AR: 


0.2000 


0.2005831 




-0.6200 


-0.6207655 




0.1520 


0.1527565 




-0.3016 


-0.3020376 


MA: 


1 .0000 


1.0001506 



We next consider a second reference model with MA order t = 2 and AR order s = 
having transfer function 

0.5 - 0 40 Jr -1 + 0.89 r~ 3 
1 - 0.20 z _1 - 0.25 z _i + 0.05 z~ 3 

The true reference model parameters and ARMA model parameter estimates are shown 
to be 



AR: 



MA: 



ACTUAL 


ESTIMATED 


0.2000 


0.1993060 


0.2500 


0.2496567 


-0.0500 


-0.0491961 


0.5000 


0.5002602 


-0.4000 


-0.3997071 


0.8900 


0.SS94749 



A third example with MA order / = 2 and AR order s = 4 having transfer function 



m 



1 + Q.2 z~’ — 0.99 z~ 2 

1 — 0.2 z -1 4- 0.62 z -2 — 0.152 z -3 4- 0.3016 z -4 



( 2 . 66 ) 
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was considered Tor which we obtained the following actual and estimated reference 
model parameter values 





ACTUAL 


ESTIMATED 


AR: 


0.2000 


0.201 1S05 




-0.6200 


-0.6223S03 




0.1520 


0.1534197 




-0.3016 


-0.3036S23 


MA: 


1.0000 


0.999763S 




0.2000 


0.199S342 




-0.9900 


-0.9SS6S52 



We consider as a final example the reference model of AR order and \1A order s = 3 and 
t — 3. respectively, with specific transfer function 



0.5 - 0.95 r~' + 1,33 z~ 2 - 0.979 z~ 3 
1 + 1.69 - -1 - 0.962 z _: + 0.2 z -3 



(2.67) 



The actual and estimated ARM A parameters are 



ACTUAL 

AR: -1.6900 
0.9620 
- 0.2000 

MA: 0.5000 
-0.9500 
1.3300 
-0.9790 



ESTIMATED 

-1.69S1325 

0.969099S 

-0.201S440 

0.4995653 

-0.9553509 

1.3346767 

-0.9S64S9S 



The above examples demonstrate the validity of the parameter estimation algo- 
rithm of [Ref. 4: pp. 619-621). Many reference models were estimated using this algo- 
rithm, including pure MA processes, for which accurate estimates were obtained. 
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D. LATTICE STRUCTURE 

In section C we developed expressions for the forward and backward prediction er- 
rors. namely, those of equations (2.31), (2.32) and (2.33). From these prediction error 
equations we can design elementary AR, MA and ARMA lattice structures or sections. 
Each elementary section satisfies the orthogonal conditions as listed in equation (2.37). 
From the prediction error recursive formulae, equations (2.31), (2.32) and (2.33), we 
construct the AR-type elementary lattice section as follows. Consider the following data 
set of order {s + 1,/) consisting of input and output data elements, 

W*) if = [ y(k) ...y(k - s) -x(k ) ... -x(k -t+ 1) -x(k - t) ] (2.68) 

W*) I 2 r = l y(k ~ 1) ...y(k -s- 1) -x(k - 1) ... -x(k - t) 0 ] (2.69) 

where If and If are the transposes of the matrices I, and I 2 defined in equations (2.55) 
and (2.56). We obtain a recursive relationship between the forward prediction errors 
v(A) of order (s + 1 j) and order ( s,t ) by substituting equations (2.52), (2.6S) and (2.69) 
in equation (2.31) such that 

v j+l,r(^) = hj+uW ^5+1,: 

= !',«.,(*) [irci+IjD^u,] (2.70) 

= v Jt XA) + d Sil {k — 1) u, 

The backward prediction error recursions are obtained in a similar manner and the 
AR-type error recursions are 

b +1 ,W = v sj(k) ~ u* djk - 1) 
l i+i,/(^) = l s,t (k) + d SJ (k - 1) 

?s + ijk) = g s /k)-u 2 d SJ (k) 

d s+\,,(k) = d 5 t {k - 1) + [ 1/3 u 3 ] v s [ (k) T - w 4 g s+u {k) 

where u, = C u\ u\ 3 and u 3 = C u$ u r 3. The AR-type elementar\ T lattice inverse section 
based on these error recursions is shown in Figure 5. 
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We design the MA-type elementary lattice inverse section in a similar manner using the 
following representation of the data set of order (sj + 1) 

h w+J ^)l3 r =[j(Ajy(A-l)...j(A-5) — .v(A') ... —x(k - It 1) -.t(A - i) ] ^ 

h w+1 (A) 1 1 = [ r(A - 1 ) ...y{k - s) 0 -x{k - 1 ) ... -x(k - i - 1 ) ] 

and by substituting equations (2.5S) and (2.72) into equation (2.31). we obtain the for- 
ward prediction error recursion as the MA order is increased by one, namely. 

£ + i (k) = vl,{k) + nUjk-\) 

l w+iM = 1 j,/ ( k ) - n \ SsJ k - 1) 

The backward prediction error relationships are obtained in a similar manner and are 
given by 

d s,t+ l(A) d s,i( k ) n 2 Ss,i 

Ss,t+ \( k ) = Ss,,( k - 1 ) - [ ”3 «3 ] Y sA k ) T ~ "A d s,i+ \( k ) 74 ^ 

The MA-type elementary lattice inverse section based on these error recursions is shown 
in Figure 6. 
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X 




We now construct the ARMA elementary lattice section from the AR and MA 
prediction error recursions. Assuming that the prediction errors are known for a given 
model order (s.t). the (s-r l.t) prediction errors can be calculated. These prediction errors 
of order (s + l.t) are then updated as the MA order increases by one resulting in a pre- 
diction error of order (s— l.t+ 1 ). We consider the forward prediction error for x as the 
AR order is increased by one. specifically 

A+i ./ (k ) = v sA k ) ~ u* li 's.t (k ~ *) (~75) 

Now. (A) becomes the current value of the forward prediction error for .y and when 
we calculate the (s+ l.t + 1 ) forward prediction error we have from equation (2.73) 

v s+u + 1 ( k ) = b+i (*) + n \ &+u (k ~ ! ) (2J °) 

Equation (2.76) can be expressed in terms of the (s.t) forward prediction errors of x by 
making appropriate substitutions for v;_ u (k) and g,_ u ,(k — 1). That is, we substitute 
qk,,, (k) and g s . u ,{k — 1) of equation (2.71) in equation (2.76) to obtain the (s+ l.t+ 1) 
forward prediction errors for .r, namely, 

W, (k) = (A) - ut d„ (A - 1) + g SJ (k - 1 )-u 2 d 5J (A - 1) ] (2.77) 
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Grouping the terms we obtain 



‘f+u +1 </(*)- + “ 2 ) l) + «f&.r(*- 1) (2.78) 

The forward prediction error recursion for y of order (s+ l,t + 1) is obtained in a similar 
manner. We begin with the (s+ l,t) order update of the prediction error and after it is 
computed, update the VIA order. Specifically, we have 

+ (2-79) 

and from equation (2.73) 

Wi (*) = &1, (*) -<&+..,(* - >) (2-SO) 

Substituting (2.71) for v>_ lr (A - ) and g s . u ,(k - 1) in equation (2. SO) then grouping terms 
we obtain the (s+ l,t+ 1) forward prediction error recursion for y. 

Wi (*) = <■ (*) + (“1 + «i « 2 ) d St , (k - 1) - ,v\ g SiI (k-l) (2.S 1 ) 

The (s+ l.t+ 1) backward prediction errors for .v and y are derived in a similar manner 
and are given by, 

Ss+i,t+\ ( k ) = Ssj ( k “ 0 + ( n 3 + "4 "3 ) (k) - («£ + th U3) l ir ( k ) 

<Wi (k) = d„ (k-l) -of vfj (A) + ttj x* t (A) 

The ARV1A elementary lattice inverse section is shown in Figure 7 where the coefficients 
are related to the prediction error recursions by the following 

ifj 1 = (zq + n* it 2 ), u'2 = n*. W3 = (zq + n\ u 2 ). iff = /q (2. S3) 

M’S = + >h It 3 ), = ("3 + «4 «3)f 1V 7 = u 3< W S = ir 3 (2.84) 






i 



X X 




Figure 7. ARM A elementary lattice section 

We see from Figure 7 that each elementary ARM A lattice inverse section contains eight 
coefTicients. 

From the AR, MA and ARMA elementary lattice inverse sections, we can obtain 
synthesis lattice structures. These structures provide a means of working with lattice re- 
alizations as linear filters. The resulting AR, MA and ARMA elementary synthesis 
lattice filters are shown in Figures S and 9 respectively. 
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Figure 8.. Top: AR elementary lattice filter. Bottom: MA elementary lattice section. 



24 



x X 




Figure 9. ARMA elementary lattice filter 

Summarizing, in this chapter we have reviewed the Mullis- Roberts criterion, intro- 
duced the ARMA parameter estimation as a generalized Mullis-Roberts criterion and 
obtained analysis and synthesis forms of lattice structures. We notice that each ARMA 
elementary lattice section consists of eight reflection parameters and the calculation of 
these parameters requires the autocorrelation and crosscorrelation information as ob- 
tained from the input output data of the reference model. Also, we obtained a set of 
equations relating the final model estimation parameters and the prediction error model 
parameters which inturn are obtained from the eight lattice parameters. 



III. ADAPTIVE LATTICE ALGORITHM 



A. LEAST MEAN SQUARE ALGORITHM 

The study and design of adaptive filters is known to be a very important part of 
statistical signal processing. Many adaptive algorithms have been developed to support 
the application of adaptive filtering in communications and control [Ref. 12). An adap- 
tive filter is characterized by the ability of its filter coefficients to adjust (self-optimize) 
automatically and yield an optimum filter design. Two processes occur within an adap- 
tive filter, namely, the adaptation and the filtering processes. During the filtering process 
a desired signal is applied to an adaptive algorithm as a reference for adjusting the filter 
coefficients. Figure 10 shows a block diagram of the adaptive modeling process. Refer- 
ring to Figure 10. let j (A) be the output of the filter at time k. By comparing the output 
with the desired signal d{k) , an error signal c(k) is generated. The adaptive algorithm 
of the filter uses this error signal to generate corrections which are applied to the filter 
coefficients such that an optimum solution is obtained. An optimization technique called 
the method of steepest descent provides an approach to solving this problem. The pro- 
cedure is as follows: 

1. Assign initial values to all filter coefficients. 

2. Using these initial values, compute the gradient vector, whose individual elements 
equal the first derivatives of the mean-squared error with respect to the filter coef- 
ficients. 

3. Compute values for the filter coefficients by changing the initial values in the di- 
rection opposite that of the gradient. 

4. Return to step 2 and repeat the procedure. 

There is. however, a limitation to this procedure. The steepest descent algorithm requires 
exact measurements of the gradient vector at each iteration which, in practice, is not 
possible. Therefore, the gradient vector must be estimated and consequently, errors are 
introduced. An algorithm is required which computes the gradient from the available 
data. The least mean square (LMS) algorithm, developed by Widrow and Hoffi is widely 
used and is very convenient to implement in real time hardware [Ref. 13: pp. 96-104). 
Let y(A) be the output of the filter and d(k) the desired signal at time k as shown in 
Figure 10. We compute the error by taking the difference between these two signals, 
namely, 



26 




Figure 10. Adaptive modeling block, diagram 



c(A-) = d(k) - y(k) (3.1) 

The value of the mean-squared error is the expected value of the error squared, 
£[[ e 2 (k) ]] and the gradient vector, V (k) , is the first derivative of the mean-squared er- 
ror. The gradient vector is given by 

where \\(k) is the time dependent filter coefficient vector. The recursion for the filter 
coefficient which changes the old value in the direction opposite to that of the gradient 
is then given by, 

„(*)_„(*_ 1) + ± - VIA)] 

e(k) 

where \v(k) is the filter coefficient vector estimate at the k ,h iteration, \\ (k — 1) is the past 
filter coefficient vector estimate, n is the convergence (gain) constant, e(k) is the error 
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tation of this algorithm proceeds as follows: 

1. Assign initial values to the filter coefficients. 

2. Compute the value for the error signal e(k). 

3. Calculate the updated estimate of the filter coefficients using the instantaneous 



4. Increment the time index by one and return to step 1. 

Convergence properties of the LMS algorithm are well documented within the literature. 
The choice of a gain constant ^ is arbitrary however, theoretical bounds have been de- 
rived for n . given by [Ref. 14: pp. 101-106], 



where / max is the maximum eigenvalue of the input autocorrelation matrix, R„. and 
where Tr [ R„ ] is the trace of the matrix R„. 

B. DERIVATION OF THE ADAPTIVE LATTICE ALGORITHM 

The adaptive lattice algorithm developed in this thesis uses concepts of the LMS 
algorithm discussed in section A and applies them to the ARM A digital lattice filter 
proposed in Chapter II. Consider the ARMA digital lattice filter of Figure 11, which 
consists of two cascaded elementary lattice sections. The filter coefficients (weights) are 
defined such that iv” represents the i' ; ' lattice coefficient at stage m of the lattice structure. 
In this figure we have a two stage lattice and there are eight coefficients per elementary 
lattice section. The output, y(k), of the lattice filter can be determined from 



gradient. 




(3.4) 



y(k) = e j{( k ) + w 4 e b 0 ( k - i) - w l %( k - *) 



(3.5) 



Forward errors at a given stage m of the lattice filter are defined as, 




(3.6) 



and the backward errors for any given stage m are, 




(3.7) 
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Figure 11. Two stage ARMA lattice digital filter. 
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where, in Figure 11. m = 1.2 and e* a = x(k) and e; = y{k). The terminal condition is 

ej m (k) = b 0 e' fm (k): m = 2 . To begin with, let b a equal unity. The initial conditions are 

e o 0 {k — 1) = 0 and e; o (k — 1) = 0. As with the LMS algorithm, we form an error between 
a desired signal </(A) and the output signal y(k) such that, 

e(k) = d(k)-?(k) (3.8) 

The instantaneous gradient according to eq (3.2) is then, 

V(A’) = 2 e(k) j^j- [ d{k) -p(k) ] (3.9) 

Since the desired signal, </(A). is not a function of the filter coefficients, equation (3.9) 

reduces to. 



V(A) - 2 c(A) [ ->“(*)] (3.10) 

where the quantitv - - r , • T — v(A) "1 is referred to as the gradient estimator. This era- 
' c\\(k) L 

dient estimator must be computed for each filter coefficient within the lattice structure. 
The filter coefficients are then updated using the respective gradient estimators. That is. 
we need to compute. 

V(A) = ^— y(k) for/= 1.2, ...,S and j= 1,2 M (3.11) 

3 ivj 

where M is the number of stages in the ARMA lattice filter. From equation (3.5), the 
gradient estimates are given by 



ern 

dwj 



Hw Swl 1 d%ik-l) cu.-] 

+ ~rj % ( k ~ ! ) + lv 4 — rr % (* ” ! ) 



>■» J 
CM; 



C \V; 



dwf 



CM; 



— US 



, *<<*-■> 



(3.12) 



cwf 



Let i jj(k) represent the partial derivatives of the output j)(A) with respect to the filter co- 
efficients and 0(A) represent the partial derivatives of the errors with respect to the filter 
coefficients. Using this representation we can re-write equation (3.12) as follows, 
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where 5H, <5j;, are kronecker deltas whose value is one if and only if 
i = 4 and j= 1 or / = 3 and j = 1 respectively. We compute (k) , by obtaining re- 
cursive relations which calculate the partial derivative of the forward and backward er- 
rors with respect to the filter coefficients. These relations are obtained from equations 
(3.6) and (3.7) by taking the partial derivatives with respect to the filter coefficients, 
namely, 



These recursive relations possess a lattice structure similar to that of Figure 11. with 
delta components injected at the summation nodes. They may also be simplified by ex- 
amining individual terms. Consider the general equation for the forward error in .v. re- 
peated here for continuity. 





(3.15) 



The partial derivative with respect to each filter coefficient is expressed as, 



Kw 8e ) 
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Since the partial derivative is taken with respect to the current filter coefficient wj at time 
A, the partial derivatives involving delay terms i.e., ( k — 1). are set to zero. This result 
follows from the realistic assumption that e x bm (k — 1) is a function of w’ ( k — 1) but not 
of ivf (k) . Also note that w> (A) is a function of w; (A — 1) but not vice versa. With these 
simplifications we reduce the equations of (3.14) to, 



4>/ m u (k) = </£_, ij (k) + b™{ (A - 1) - <5?/ 4 m _ t (k- 1) 

u w = c, u (a> + 1); < (a ■ - 1) ■ - < (* - 1) 

< u (A) = 77 (A) + w s C-, u 7') - 77 7L (7 - < <t>L, ij (A) 



(3.17) 






(A) = - <7/ (A) - C, V (A) + 5?/ (A) + vv-r C, ly (A) 



and the gradient estimator is. 



7 j (A) = 0/, / j (A) + 









(3 IS) 



Although these arc valid recursive relations, they are difficult to implement in a lattice 
algorithm. The ultimate goal is the requirement to easily compute i p,j{k) from the 
available data. From eq (3.18) it is evident that ^ i} ( k ) depends on 4>} VJ ( k ) which in turn 
requires knowledge of ( k ) and (A) but not of <f>l m , 7 (A) or <f>i miJ (A). Therefore the 

three equations necessary to compute the gradient estimator are. 

hj (A) - ^ /; (A) + 6\\ ej (A - 1) - 4 0 (A - 1) 

^iy(A) = C , iy (A) + #/<_, ^ " 1) - (A - 1) (3.19) 

€o(A) = 0L li7 (A) + ^; +1); <(A- \)-5f^ )j 4 m {k- 1) 



These equations are dependent on the filter coefficients wj, / = 1 .2.3.4 and 

j = 1,2 M , thereby reducing by one-half the number of computations required for 

4>' mlJ (A) and (A). A recursive relation is desired for \p,j(k) which does not involve 
delta functions. Consider the four stage lattice filter with terminal condition b 0 equal to 
unity such that (A) = (A) . The procedure for computing ^, V (A) using the 

equations of (3.19) is as follows: 

1. Calculate 0/ <,-(A) with m equal to one and letting i and j range from one to four. 
Repeat for m = 2,3,4. 

2. Using the terminal condition and expression for <^ m , v (A), calculate </>) ](; (A). 

This requires solving 88 equations, however, the result is a very simple recursive formula 
for i jj tj (A) , namely, 
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/ = 1 . 3 



(3.20) 



= (*- J ) 

^ ij {k) = e x b) Jk- 1) /= 2, 4 

The lattice coelTicients are calculated by substituting the recursive formula of eq (3.20) 
for the gradient estimator in equation (3.3). The adaptive coefficient update equation is, 

>4(*) = »4(A-1)-MA)iM*) (3-2i) 

Since the gradient estimator, and therefore the gradient, is the same for vwf , / = 1,3 and 
w’, i = 2,4 , it follows that w{ = wj and \v’ 2 = wj. The number of filter coefficients re- 
quired to update the lattice filter is reduced to two, i.e., vv{ and w' 2 . Furthermore, from 
the symmetry of the lattice structure, the following equalities between filter coefficients 
are assumed, 



u-i = 
w i = 



ll 5 

*7 

4 



= 

vr( = vr£ 



(3.22) 



Incorporating these equalities with those derived by the gradient estimator produces the 
elementary ARMA lattice section of Figure 12 where, 

u 2 = u 4 = M'S = 11 6 = rj 
j j j j l 

= 1 '3 = u 7 = = k j 



To prove that these coefficient reductions are valid, a computer generated solution using 
the Fortran program of Appendix D was compared to hand analysis of a second order 
transfer function and lattice filter. The output of the ARMA digital lattice filter was first 
put into difference equation form and then compared to the known transfer function. 
From this comparison, lattice coefficients were computed. Details of this analysis are as 
follows. Consider a two stage ARMA digital lattice filter comprised of the reduced ele- 
mentary section shown in Figure 13 and a transfer function of the form. 




b 0 + 6, 2 + bj z 

1 + r?i z + O.J 2 



(3.24) 
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Figure 12. Simplified elementary ARMA lattice section. 



The output of the lattice filter can be written in difference equation form by earn ing out 
the following steps: (i) start with the output of the lattice filter equation (3.5), (ii) sub- 
stitute expressions for the forward and backward errors, equations (3.6) and (3.7). re- 
spectively. into equation (3.5) and (iii) earn out the algebra. A detailed derivation of this 
difference equation is given in Appendix E. The difference equation in its final form is 
given by, 



y(k) — x(k) + 2( vcj + ivj wf + 



— 2 ( w} + 



i 



uy 



2 

vv*2 ■+■ 



ivy uy )x(k — 1)4-2 ivy x{k — 2) 
vv,' uf M* ~ 1) ~ 2 w]v{k - 2) 



(3.25) 



which can be written in the transfer function form as, 
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(3.26) 



ii / -w 1. 1 2. 1 2\ — 1 . ^ 2 — 2 

1 4- 2( w 2 4- u*j ;v*j 4- w 2 w 2 ) z 4-2 \v 2 z 

T~T7f i ! i 2 ! 1 2 x -1 . ~ 2 -2 

1 4- 2( iv, + ir 2 vi '2 + iv, iv, ) z +2 iv, z 



Comparing the lattice filter transfer function, 7//z) with the known filter transfer func- 
tion // ; (z), produces the following relationships between filter coefficients, 



l -)/->, 12 12s 

£>, = 2( iv 2 + iv, iv, + rv 2 vv 2 ) 
a, = 2( iv, + iv 2 iv 2 + iv, iv, ) 
b 2 = 2 iv*2 
a 2 — 2 iv. 



(3.27) 



Solving for ivf and iv; in terms of the known transfer function coefficients b 2 and a 2 and 
then substituting these results into the expressions for b x and respectively, yields the 
following. 



b\ = a 2 w \ + ( 2 + b 2 ) ivs 
n, = ( 2 4- a 2 ) iv, 1 + b 2 n '2 



or in matrix form, 

i 

IV, 
i 

IV, 



" V 




a 2 


2 + b 2 




. c7 l _ 




_2 + a 2 


b 2 





(3.2S) 



(3.29) 



and solving for the lattice coefficients, we have 



1 


b 2 


— (2 + b 2 ) " 


~b,~ 


a 2 b 2 - (2 4- b 2 ) (2 4- a 2 ) 


_ - (2 + a 2 ) 


a 2 


. a \ . 



and 



IV, 




« 2 
1 2 



(3.31) 



Now that a method of converting between lattice and transfer function coefficients for 
a second order system has been established, we consider the specific transfer function 
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H/z) = 



1 - O.S z _1 + 1.78 z -2 



(3.32) 



1 - 0.S9 z -1 + 0.25 z -2 



where 



* 0 -1.0 b x = —0.80 * 2 =1.78 

= —0.89 a 2 = 0.25 



From equations (3.30) and (3.31) the lattice coefficients are calculated as, 



w 2 = 0.125, w 2 = 0.890, w} = -0.240719, w\ = -0.195719 



Values for the steady-state lattice coefficients were computed using the Fortran program 
in Appendix D and are shown below. Convergence aspects of both the lattice coeffi- 
cients and output error are shown in Figure 13. 



these results confirm the validity of the derived adaptive lattice algorithm and the design 
of a new elementary lattice section shown in Figure 12. 

The current adaptive lattice algorithm assumes that the terminal condition is unity. 
This is generally not the case in practice. We now extend this adaptive algorithm to the 
more general case where the terminal condition is an arbitrary constant. The recursive 
relation which updates * 0 is similar to those which update the other lattice filter coeffi- 
cients. The update equation for b 0 is given by 



the desired signal is not dependent on b 0 . The gradient estimator for b 0 is written 
as, 



Since the partial derivative is taken with respect to b 0 at time k, this reduces to, 



w 2 = 0.1249S2, w\ =0.890003, w{ =-0.240710, w 2 ' =-0.195711 




(3.33) 



The gradient estimator 
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cv(*) de )-ik) 

cb 0 £b 0 



Similarly, 

e L (*) = i-i M + vy 4 ,+1 e b m ( k - !) - » 3 m+1 4 m ( k - 1 ) (3-36) 



and the partial derivative with respect to b 0 is 




cb 0 db 0 



The terminal condition is. 






[k) = b { 



er..(k) 



0 fc / w 



(3.37) 



(3-3S) 



Taking the partial derivative of equation (3.38) with respect to b 0 yields e? w (A), and the 
recursive equation to update b 0 ( k ) becomes, 

b 0 (k) = b 0 (k-l)-ne(k)e? M (k) (3.39) 



The gradient estimators for the lattice filter coefficients are scaled by the arbitrary con- 
stant b 0 . since at the terminal condition (A) = b 0 (A) and b 0 is propagated 
through the calculations. The gradient estimators become, 



'Pij ( k ) — ~ bo j (A — 1) / — 1 ,3 
t iJ (k) = b 0 c x bj i {k- I) i = 2.4 



To test this more general adaptive lattice algorithm the output of a known transfer 
function with b 0 equal to 0.5 was compared to the output of the ARM A digital lattice 
filter. The second order transfer function used was, 



Hj{2) 



0.5- 0.4 z~' + 0.S9 2~ 2 
1 — 0.S9 z -1 + 0.25 z~ 2 



(3.41) 



The computer generated steady state lattice coefficients are given below and convergence 
aspects shown in Figure 14. 

b 0 = 0.499946, w' 2 = -0.120428. vv 2 2 = 0.593237, w] = -0.447423, wf = 0. 166320 
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Figure 14. Top: Lattice coefficients. Bottom: Output error. 
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In order to maintain the same adaptive time constant and misadjustment at each 
stage in the lattice, the convergence constant is normalized by the power level at each 
stage [Ref. 15]. Therefore, we can write equations (3.21) and (3.39) as, 

4 (*) = M (* - 1) - -T77 - m t tJ ( A ) 

Oj (A) 

1 (3-42) 

b Q (k) = b 0 (k-l)--f— e{k)e} M (k) 

where p is the convergence constant and aj (k) and y 2 (A) are estimates of the power at 
the j' h stage for vv' and b 0 respectively and computed as follows: 

°) (A) = P o) (A - 1) + (1 - p) tfj(k) 

y 2 (A) = p y 2 (A - 1) + (1 - p) [ (A) ] 2 (3 ' 4j) 

Writing equation (3.42) using the notation adopted for the reduced elementary ARMA 
lattice section we obtain 



>)(*)“'>(*- 1 ) T~7~ 

Cj (A) 

kj (A) = kj (A — 1) —— 

o] (A) 

b 0 (k) = b 0 (k-l)--r— 

r (A) 



e (A) e bj_ x (A ~ 1) 

(A-i) 



e(k) e/ M (A) 



(3.44) 



In the above equations p is a weighting parameter, 0 < p < 1, which distributes the 
amount of weight given the past power level or current sample. Normalized convergence 
constants are used in all examples of this thesis. The adaptive lattice algorithm is sum- 
marized in Table 1. 

In summary, we have derived an adaptive algorithm based on the LMS theory of 
adaptive coefficient computation. This new adaptive algorithm easily updates the lattice 
coefficients by using available data. The original requirement to update eight coefficients 
of an elementary ARMA lattice section was reduced to updating only two coefficients 
and still being able to describe the lattice. The algorithm is general in that it applies to 
systems whose terminal condition is an arbitrary constant. The validity of this algorithm 
was demonstrated through comparisons between hand analysis and computer simu- 
lation. In the next chapter, we further demonstrate the convergence of this algorithm. 
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Table 1. SUMMARY OF ADAPTIVE LATTICE ALGORITHM 

With given initial conditions, e x h = x{k), ef g = y{k), ef, 0 (k - 1) = c x bQ (k - 
lattice coefficients zero, 

Step 1: for k= l,m compute 



£<*> - ‘L> <*) + < (* - 1) 

i (*) = (*) + »" + ‘ <(*- d - < + ’ < (* - 1 ) 



with output e>- g (k) 

Step 2: for k= l,m compute 



Step 3: 



J m— 1 






(*) 

(*) 



Update coefficients 



r j (k) = fj (k — 1) — 

A y (/:) = (A - 1)- 
* 0 (A-) = ^ 0 (* - 1) - 




e{k)e x bj i {k — 1) 

4 *)^., (*- 1 ) 




4 *) % (*) 



Step 4: Repeat for next iteration i.e. return to step 1. 



= 0 and all 
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IV. EXPERIMENTAL RESULTS 



The adaptive lattice algorithm derived in Chapter III is now computer simulated to 
study its convergence performance. The system identification mode of adaptive filtering 
is considered for this purpose. Figure 15 shows a general system identification config- 
uration. The systems considered are time-invariant and linear. Notice that we apply the 
same input, white noise in general, to both the reference system and the adaptive lattice 
filter which is modeling the system. The criterion in this configuration is to minimize the 
mean-squared error between the system and filter outputs. Thus, in this context, the 
adaptive algorithm continuously updates the lattice filter parameters in order to mini- 
mize the mean-squared error. 

The adaptive algorithm is realized as summarized in Table 1. As we mentioned in 
Chapter III. the two important parameters of the algorithm are the adaptation constant 
p and the weighting constant p . In what follows, we shall consider convergence studies 
of both second and third order reference systems (fixed filter transfer functions). Con- 
sider the following reference system with transfer function, 



Uj{z) 



1 + 0,2 z -1 — 0.35 z~ 2 
1 - 1.4 z -1 + 0.S5 z~ 2 



This system has complex poles and simple zeros located at z = (0.7 ±j0.6) and 
z = 0.5, —0.7, respectively. Using a convergence constant p = 0.01 and power level 
weighting factor p = 0.45, the adaptive ARMA digital lattice filter which models the 
above system has the following steady-state lattice parameters. 



terminal condition b 0 = 0.999202 

lattice coefficients r, 1 = 0.3525S0 
r, 2 = -0.174355 
kj = — 0.44S228 
k?= 0.425060 

Convergence properties of the lattice coefficients and error are shown in Figure 16. The 
mean-squared error was minimized after approximately 1700 iterations at which time the 
lattice coefficients reached their steady-state values. When the value of the convergence 
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Figure 15. Block diagram of system ^identification modeling 



constant p was modestly increased, convergence degraded rapidly. Also, when the 
weighting factor p was increased, convergence deteriorated quickly. From this, we con- 
clude that the convergence constant is the more sensitive input parameter. 

Let us consider another second order dynamic system with transfer function, 



//,(-) 



0.5 - 0.2 z -1 - 0.445 r~ 2 
1 -z" 1 + 0.94 z -2 



This system has complex poles at z = (0.5 ±jQ.S) and complex zeros located at 
z = (0.2 ±j0.7). Using a convergence constant p = 0.005 and power level weighting 
factor p = 0.97, the adaptive ARMA digital lattice filter which models this system has 
steady-state lattice parameters as follows, 



terminal condition b 0 = 0.500027 

lattice coefficients r] = 0.104654 
r] = 0.296658 
k\= -0.429232 
k\ = 0.627718 



43 




44 



Convergence properties of this adaptive filter are shown in Figure 17. In this example, 
convergence was obtained after approximately 2500 iterations. Again, by changing the 
values of p and p slightly, covergence deteriorated with the convergence constant p being 
the more sensitive parameter. 

Next we consider a third order reference system with known transfer function, 



///;) = 



0.5 - 0.95 z~‘ + 1.33 z~ 2 - 0.979 z ~ 3 
1 - 1.69 z -1 +0.962 z -2 - 0.2 z“ 3 



The adaptive ARMA digital lattice filter which describes this system has the following 
steady-state lattice parameters, 

terminal condition b 0 = 0.499970 

lattice coefficients rj = —0.328447 
r] = 0.399472 
r\ = -0.652706 
k\ = — O.S2173S 
k] = -0.091111 
A?= -0.133333 



These parameters were obtained using a convergence constant p = 0.015 and power level 
weighting factor p = 0.9. Convergence properties of this adaptive filter are shown in 
Figure 18. Steady-state values for the lattice coefficients were obtained after approxi- 
mately 7100 iterations. It is reasonable to assume that a third order system will converge 
more slowly than a second order system. The number of iterations required for this third 
order system to converge is consistant with convergence rates of other adaptive algo- 
rithms which model third order systems [Ref. 16]. The input parameter p was again 
found to be the more sensitive parameter. 

In all the previous examples, the values of p and p may or may not be optimum 
values. That is, an exhaustive search of all combinations of p and p was not performed 
to demonstrate convergence of the algorithm. Nevertheless, a number of different ways 
of realizing the value of the convergence constant p have been reported in the literature. 
In one method, Mikhael et. al. [Ref. 17] have obtained a variable p by using a self opti- 
mizing technique. In this method, p is calculated from the input data as an iteration 
process and is individually determined for each filter parameter. In another method p is 
chosen by using a variable step LMS technique [Ref. 18], where the range of p is 
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Figure 18. Third order ARMA lattice filter, terminal condition of 0.5. 
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specified by p max and p mm which are within the bounds described b\ equation (3.4) of 
Chapter III. These techniques of choosing p during the adaptation process have been 
shown to improve filter convergence. They, however, require additional computations 
to achieve this faster convergence. When a combination of the parameters p and /> was 
obtained which yielded convergence, these values were chosen for examples. Besides the 
examples reported, simulation studies have been carried out for several other cases. In 
all cases, however, definite convergence of the algorithm has been observed. 

In summary, we have demonstrated through computer simulation that the derived 
adaptive lattice algorithm is suited for system identification modeling. Furthermore, we 
have shown that there is flexibility in choosing the values of the convergence constant 
p and weighting factor p . Some techniques for selecting (computing) the value of the 
convergence constant have been introduced. These methods improve convergence at the 
cost of additional computations. 

A. CONCLUSIONS 

In this thesis we have demonstrated that the ARMA parameter estimation algo- 
rithm proposed in [Ref. 4: pp. 619-621] is a valid method for obtaining approximations 
to reference models. Furthermore, the criterion used to derive the algorithm is a gener- 
alized form of the Mullis-Roberts criterion for least squares modeling. The AR and MA 
parameters of the ARMA model can be updated independently as their respective orders 
increase by one. From the recursive prediction error formulas, an ARMA digital lattice 
filter was designed with arbitray AR and MA orders. 

For the ARMA digital lattice filter, we derived an adaptive lattice algorithm. This 
algorithm was based on the least mean square method of optimizing coefficients. The 
derived adaptive lattice algorithm can easily compute the values of the lattice coefficients 
from available data. The algorithm simplified the number of coefiicents required to be 
updated from eight coefficients per elementary lattice section to only two such that the 
filter can be completely described. This savings in computational effort makes the al- 
gorithm attractive for identification of unknown systems since many systems require an 
ARMA model for parsimonious modeling. 

Convergence of the adaptive lattice algorithm was demonstrated with several exam- 
ples in Chapter III. The number of iterations required before convergence varied greatly 
between second and third order models as well as within second order models. Optimum 
convergence rates were not sought after as much as proving the convergence of the al- 
gorithm. Rapid convergence rates were demonstrated in Chapter II for a second order 
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system upon completion of an extensive search for the optimun values of the conver- 
gence constant and power level weighting factor. 

Although a method of converting between direct form and lattice realizations for 
second order ARM A filters was developed, a general algorithm to perform this trans- 
formation given the ARMA lattice filter design was not obtained. When solving for a 
transformation between filter realizations of third order the solution is hindered by 
nonlinearities. 

The objectives of the thesis were sucsessfully accomplished. Some suggestions for 
future work include the following: (i) extensive theoretical analysis for determining op- 
timum values for n , (ii) derivation of a generalized algorithm which converts any given 
ARMA transfer function into a set of lattice parameters, (iii) development of theoretical 
convergence models for the ARMA adaptive lattice algorithm and analysis of these 
models and (iv) application of the adaptive lattice filter, both analysis and synthesis 
forms, in modeling such practical signals as speech. ARMA lattice filter modeling has 
considerable application potential because of its very accurate modeling of nearly any 
signal or system of interest. 
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APPENDIX A. MAIN PROGRAM TO ESTIMATE ARMA PARAMETERS 



c 
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c 



THIS PROGRAM COMPUTES THE ARMA PARAMETERS AS THE AR OR MA ORDER 
OF AN ARMA MODEL INCREASES BY ONE. IT USES THE ARMA PARAMETER 
ESTIMATION ALGORITHM PROPOSED BY MIYANAGA, NAGAI AND MIKI. 

VARIABLE DEFINITIONS 



VN - INPUT VECTOR CONTAINING DATA GENERATED AS WHITE GAUSSIAN 
NOISE 

Y - OUPUT VECTOR OF DIFFERENCE EQUATION WITH VN AS INPUT 
RX - AUTOCORRELATION DATA OF INPUT VN 

RY - AUTOCORRELATION DATA OF OUTPUT Y 

RXY - CROSSCORRELATION DATA OF INPUT AND OUTPUT 

RYX - CROSSCORRELATION DATA OF OUTPUT AND INPUT 

NDATA - NUMBER OF INPUT DATA POINTS 

KDATA - NUMBER OF INITIAL DATA POINTS TO DISREGARD 
XI - INPUT DATA VECTOR AFTER DISCARDING KDATA POINTS 
Y1 - OUTPUT DATA VECTOR AFTER DISCARDING KDATA POINTS 
A - VECTOR CONTAINING CALCULATED AR PARAMETERS 
B - VECTOR CONTAINING CALCULATED MA PARAMETERS 
TXA - VECTOR CONTAINING COEFFICIENTS FOR FORWARD PREDICTION OF 
INPUT 

TYA - VECTOR CONTAINING COEFFICIENTS FOR FORWARD PREDICTION OF 
OUPUT 

TXB - VECTOR CONTAINING COEFFICIENTS FOR FORWARD PREDICTION OF 
INPUT 

TYB - VECTOR CONTAINING COEFFICIENTS FOR FORWARD PREDICTION OF 
OUTPUT 

GA - VECTOR CONTAINING COEFFICIENTS FOR BACKWARD PREDICTION OF 
INPUT 

GB - VECTOR CONTAINING COEFFICIENTS FOR BACKWARD PREDICTION OF 
INPUT 

ZTA - VECTOR CONTAINING AR COEFFICIENTS FOR BACKWARD PREDICTION 
OF OUTPUT VALUES 

ZTB - VECTOR CONTAINING MA COEFFICIENTS FOR BACKWARD PREDICTION 
OF OUTPUT VALUES. 

NI - LENGTH OF DATA VECTOR XI 

NK - LENGTH OF DATA VECTOR XI MINUS ONE USED TO START 

CORRELATION COMPUTATIONS. 

NS - DESIRED AR ORDER OF ARMA MODEL. 

NT - DESIRED MA ORDER OF ARMA MODEL 

KS - CURRENT AR ORDER OF UPDATE 

KT - CURRENT MA ORDER OF UPDATE 

VX - EXPECTED VALUE OF PREDICTION ERROR FOR INPUT SQUARED. 

VY - EXPECTED VALUE OF PREDICTION ERROR FOR OUTPUT SQUARED. 

VXY - EXPECTED VALUE OF PRODUCT OF PREDICTION ERROR FOR INPUT 
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o n n 



c 

c 

c 

c 

c 

c 

c 



AND OUTPUT. 

VG - EXPECTED VALUE OF BACKWARD PREDICTION ERROR OF INPUT 
SQUARED. 

VZ - EXPECTED VALUE OF BACKWARD PREDICTION ERROR OF OUTPUT 
SQUARED. 

VGZ - EXPECTED VALUE OF PRODUCT BETWEEN BACKWARD PREDICTION 
ERRORS OF INPUT AND OUTPUT. 

DIMENSION VN( -5: 2006) ,Y(-5: 2006) ,RX(0: 2000) ,RY(0: 2000) ,RXY( 0: 2000) 
DIMENSION RYX( 0: 2000) ,X1( 2000) ,Y1( 2000) ,X( 12) ,A( 0: 21),B(0: 21) 
DIMENSION TXA( 0: 21) ,TYA(0: 21) ,TXB(0: 21) ,TYB(0: 21) ,GA(0: 21) 
DIMENSION GB( 0: 21) ,ZTA(0: 21) ,ZTB(0: 21) 

INPUT DATA INFORMATION 



1 



2 

C 

C 

C 

C 

C 



10 



15 



20 



WRITE (6,1) 

FORMAT (/’ ENTER THE NUMBER OF DATA POINTS: ') 

READ (6,*) NDATA 
WRITE (6,2) 

FORMAT (/' ENTER NUMBER OF INITIAL DATA POINTS TO DISREGARD: ') 
READ (6,*) KDATA 



•jWr-jVvWrvWr 5W0V* VrsWoWrVcV?** VcVoWoV VnV Vr ■ 



’ ■}’>- ‘V Vc ■} V Vr Vc Vr 'k Vr * Vr ■>> Vr Vr 'V Vr Vr ■>> vV ‘/r -jV ■> V Vr ■ 



■jWrVc'WrV-V-’ 



Vr-ZcsWr 



INITIALIZE ARRAYS 

DO 10 L=-5, NDATA 

VN( L)=0 

Y(L) =0 

CONTINUE 

DO 15 L=0 , 20 

A(L)=0 

B(L)=0 

TXA(L)=0 

TXB(L)=0 

TYA(L)=0 

TYB(L)=0 

GA(L) =0 

GB( L) =0 

ZTA(L)=0 

ZTB( L)=0 

CONTINUE 

DO 20 L=0 , 1999 

RX(L)=0 

RY( L)=0 

RXY(L)=0 

RYX(L)=0 

CONTINUE 

DO 25 L=1 , 12 
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X(L)=0 

CONTINUE 



25 

C 

C 

C 

C 

c 



35 

30 

C 

C 

C 

C 

c 

c 

c 

c 

c 

c 

c 

c 

c 









■VoWr-V 



■ Vc">V V? Vc V? < V V? Vc Vr ‘V Vr V? Vr V? ■> V sV Vc V? Vr V? Vc Vr ‘V V? “ V “ wV~V < V Vr Vr Vr Vc 



GENERATE WHITE GAUSSIAN NOISE INPUT 

ISIZE=NDATA 

IX=152255 

ISORT=0 

MUL=2 

DO 30 K=1 , ISIZE 

CALL SRND( IX, X, 12 ,MUL, ISORT) 
XT=-6. 0 

DO 35 1=1,12 
XT=XT+X( I) 

VN(K)=XT 

CONTINUE 



■> V Vc Vc -jV Vr Vc :V Vr Vr Vc Vc Vc Vr Vc Vr Vr Vc ■> V Vc * V? Ve Vc Vc * “V Ve Vc * Vr Vc Vc Vc x Vr V? * Vc V? Vc Vc •>> V? Vc VoV “V "k Vf Vc VrtV Vr * Vr Vr Vr Vc V? Vr sV Vc Vr V c Vc 



COMPUTE OUTPUT OF REFERENCE MODEL FILTER AND DISREGARD SPECIFIED 
NUMBER OF DATA POINTS 

DO 40 L=1 ,NDATA 

Y(L)=VN(L)+1. 6*Y( L- 1) -0. 95*Y(L-2) 

Y(L)=VN(L)+0. 2*Y(L-1) -0. 62*Y(L-2)+0. 152*Y( L-3) -0. 3016*Y(L-4) 

Y( L)=VN( L) -0. 2*VN(L-l)+0. 62*VN( L-2) -0. 152*VN( L-3)+0. 3016*VN(L-4) 

Y( L)=VN( L)+0. 2*VN(L-1) -0. 99*VN( L-2)+0. 2*Y(L-l)-0. 62*Y(L-2)+0. 152*Y 
&(L-3)-0. 3016*Y(L-4) 

Y( L)=VN( L) - 1. 6*VN(L- 1)+1. 45*VN( L-2)+l. 2*Y(L-l)-0. 72*Y(L-2) 
Y(L)=VN(L)+0. 2*VN(L-l)-0. 35*VN( L-2)+l. 4*Y(L-l)-0. 85*Y(L-2) 

Y( L)=0. 5*VN( L) -0. 2*VN(L-l)+0. 445*VN(L-2)+Y(L-l) -0. 94*Y(L-2) 



C Y( L)=VN( L) -2. 7*VN(L-l)+3. 21*VN(L-2) -1. 595*VN(L-3)+l. 95*Y(L-1)-1. 62 
C &*Y(L-2)+0. 54*Y(L-3) 

C Y( L)=VN( L) - 1. 0*VN( L- 1)+0. 89*VN( L-2)+0. 40*Y(L-l)-0. 2121*Y( L-2) -0. 20 
C &894*Y(L-3)-l. 810373*Y( L-4) 

C Y( L)=0. 5*VN(L) -0. 95*VN(L-1)+1. 33*VN( L-2) -0. 979*VN( L-3)+l. 69*Y(L-1) 
C &-0. 962*Y( L-2)+0. 20*Y(L-3) 

C Y(L)=0. 5*VN( L) -0. 4*VN(L-l)+0. 89*VN( L-2)+l. 69*Y(L-l)-0. 962*Y( L-2)+0 

C &. 2*Y(L-3) 

C Y( L)=0. 5*VN(L) -0. 4*VN(L-l)+0. 89*VN(L-2)+0. 2*Y(L-l)+0. 25*Y(L-2)-0. 0 

C &5*Y( L-3) 

C Y(L)=0. 5*VN( L) -0. 4*VN(L-l)+0. 89*VN( L-2)+0. 89*Y(L-l)-0. 25*Y(L-2) 

C Y(L)=. 0154*VN(L)+. 0642*VN(L-l)+0. 0642*VN(L-2)+0. 0154*VN(L-3)+l. 99* 

C &Y(L-1) -1. 57*Y( L-2)+0. 4583*Y(L-3) 

C Y(L)=0. 5*VN(L)+0. 256*VN( L-l)+0. 1234*VN(L-2)+0. 0987*VN(L-3) 

40 CONTINUE 
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LJ=KDATA+1 
DO 45 L=LJ , NDATA 
LK=L-KDATA 
X1(LK)=VN(L) 

Yl( LK)=Y(L) 

CONTINUE 

WRITE (*,77) (Y1(K), K=200,1800) 
FORMAT (5( 1X,F10. 6)) 



Vc Ve 7 V Vr Vc Vc Vr V? -jV Vr Vc Vc -jV sV Vc 7 V * Vr * * Vc Vc * Vc *V -5 V Vc * Vc “V Vr 7V Vc Vc Vc ?Y 7V Ve Vr 7 V 7 V 7V -jV Vc Vr :V Vr Ve *V 7V 7V * Vc Vr 7 'c *jY Vr Vc V r Vr ~,V 



COMPUTE AUTO-CORRELATION AND CROSS -CORRELATION TERMS 

NI=NDATA-KDATA 
NK=NI- 1 

CALL CORREL (NI , 50 ,X1 , Y1 ,RX,RY,RXY,RYX,NK) 

DO 50 L=0 , 10 

WRITE (*,200) RX(L) ,RY(L) ,RXY(L) ,RYX(L) 

WRITE (9,200) RX(L) ,RY(L) ,RXY(L) ,RYX(L) 

WRITE (9,201) 

CONTINUE 
WRITE (9,201) 

FORMAT (2X,4(2X,F14. 9)) 

FORMAT (’ ') 



^Wf'sWrsWrsWoWoV’jWrjWr'sWoWrjV'sV'sWoV'sWc^WrsV’jV'jV^V'sV'sV' 



tHWc: sWrVcVf’sV’jWrVr'sWfVoWoV’sWf^WrVrVc’jWr’sV 



INPUT THE DESIRED AR AND MA ORDERS THEN DEFINE INITIAL CONDITIONS 
WRITE (6,3) 

FORMAT (/' ENTER THE DESIRED AR ORDER: ') 

READ (6,*) NS 
WRITE (6,4) 

FORMAT (/’ ENTER THE DESIRED MA ORDER: ’) 

READ (6,*) NT 



KS=0 

KT=0 

A(0)=1. 0 
VX=RX( 0) 
VY=RY( 0) 
VXY=-RYX( 0) 
VG=RX( 0) 
VZ=RY( 0) 
VGZ=-RYX( 0) 
TXB( 0)=1. 0 
TYA( 0)=1. 0 
GB( 0) =1. 0 
ZTA( 0)=1. 0 
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VrVr VrVr VrVr VrVr VrVr VrVr VrVrVr VrVrVr VrVr VrVr Vr VrVrVr' 



VrVr VrVr VrVrVr VrVr * VrVrVr * VrVr 



• VrVrVr VrVr VrVrVr Vr Vr Vr Vr 



• Vr Vr Vr Vr Vr Vr Vr V r Vr 



ESTIMATE THE ARMA PARAMETERS 



300 



301 



C 

302 



303 



C 

C 

C 

C 

C 



211 



IF (NT. EQ.0. AND.KS.LT. NS) THEN 
KS=KS+1 

CALL NEWAR(KS , KT , VX , VY , VXY , VG , VZ , VGZ , TXA , TXB , TYA , TYB , GA , GB , ZTA , ZTB 
& , RX , RY , RXY , RYX , A , B ) 

GOTO 300 
ELSE 

IF (NS. EQ. 0. AND.KT.LT. NT) THEN 
KT=KT+1 

CALL NEWMA(KS , KT , VX , VY , VXY , VG , VZ , VGZ ,TXA , TXB , TYA , TYB , GA , GB , ZTA , ZTB 
& , RX , RY , RXY , RYX , A , B ) 

GOTO 301 

ENDIF 

ENDIF 

IF (NS.NE. O.OR. NT. NE. 0) THEN 
IF (NS. GE. NT. AND. NT. NE. 0. AND. KT. LT. NT) THEN 
KS=KS+1 

CALL NEWAR( KS , KT , VX , VY , VXY , VG , VZ , VGZ , TXA , TXB , TYA , TYB , GA , GB , ZTA , ZTB 
& , RX , RY , RXY , RYX , A , B ) 

KT=KT+1 

CALL NEWMA( KS , KT , VX , VY , VXY , VG , VZ , VGZ , TXA , TXB , TYA , TYB , GA , GB , ZTA , ZTB 
& , RX , RY , RXY , RYX , A , B ) 

GOTO 302 
ELSE 

IF (KS.LT. NS) THEN 
KS=KS+1 

CALL NEWAR( KS , KT , VX , VY , VXY , VG , VZ , VGZ , TXA , TXB , TYA , TYB , GA , GB , ZTA , ZTB 
& , RX , RY , RXY , RYX , A , B ) 

GOTO 303 
ENDIF 
ENDIF 
ENDIF 



Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr 



Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr it Vr Vr Vr Vr Vr Vr Vr Vr Vr Vr • 



VrVrVrVr 



PRINT ESTIMATED ARMA PARAMETERS 

WRITE (*,211) 

WRITE (9,211) 

WRITE (*,210) (A(K) , K=1,KS) 
WRITE (9,210) ( A(K) , K=1,KS) 
WRITE (*,211) 

WRITE (9,211) 

F0RMAT(’ ’ ) 

WRITE (*,210) (B(K) , K=0,KT) 
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210 



C 

C 

C 

C 

C 



60 



70 



WRITE (9,210) ( B(K) , K=0,KT) 
FORMAT (' ' ,1X,4(2X,F13. 10)) 
STOP 
END 



•>V Vf *V VvVr 



■Vr *sV }Wc ?W«V Vo’oWoV Vc Vr VoV 1 



■ *}V ->V -j V ■> V * *>V * Vc Vc Vc Vc •jV Vr ■; V Vc ‘j V Vc Vr ■sV Vr Vc Vc V- ■> V Vc “V Vr * •> V Vr Vr ■>> ‘V Vr Vc 



SUBROUTINE TO COMPUTE CORRELATION TERMS 

SUBROUTINE CORREL(N , LAG , X , Y , RX , RY , RXY , RYX , NK1 ) 

REAL X(0: NK1) ,Y(0: NK1) ,RX(0: 2000) ,RY(0: 2000) ,RXY( 0: 2000) 
REAL RYX( 0: 2000) , SUM1 ,SUM2 , SUM3 , SUM4 
DO 70 K=0 , LAG 
NJ=N-1-K 
SUM 1=0 
SUM2=0 
SUM3=0 
SUM4=0 
ANK=NJ 

DO 60 J=0 , N J 

SUM1=SUM1+X(J+K)*X(J) 

SUM2=SUM2+Y(J+K)*Y(J) 

SUM3=SUM3+X( J+K)*Y( J) 

SUM4=SUM4+X( J)*Y( J+K) 

CONTINUE 
RX( K)=SUM 1/ANK 
RY(K)=SUM2/ ANK 
RYX(K)=SUM3/ANK 
RXY(K)=SUM4/ANK 
CONTINUE 
RETURN 
END 
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APPENDIX B. SUBROUTINE FOR MAIN PROGRAM 



c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



SUBROUTINE NEWAR( KS , KT , VX , VY , VXY , VG , VZ , VGZ , TXA , TXB , TYA , TYB , GA , GB , Z 
&TA , ZTB , RX , RY , RXY , RYX , A , B ) 

THIS SUBROUTINE COMPUTES AR PARAMETER VALUES FOR AN ARMA MODEL AS 
THE AR ORDER INCREASES BY ONE. 

VARIABLE DEFINITIONS 

KS - CURRENT AR ORDER 

KT - CURRENT MA ORDER 

RX - AUTOCORRELATION DATA OF INPUT VN 

RY - AUTOCORRELATION DATA OF OUTPUT Y 

RXY - CROSSCORRELATION DATA OF INPUT AND OUTPUT 

RXY - CROSSCORRELATION DATA OF OUTPUT AND INPUT 



c 

c 


TXA 


- VECTOR CONTAINING 
OF INPUT. 


AR 


COEFFICIENTS 


FOR 


FORWARD 


PREDICTION 


c 

c 


TXB 


- VECTOR CONTAINING 
OF INPUT 


MA 


COEFFICIENTS 


FOR 


FORWARD 


PREDICTION 


c 

c 


TYA 


- VECTOR CONTAINING 
OF OUTPUT 


AR 


COEFFICIENTS 


FOR 


FORWARD 


PREDICTION 


c 

c 


TYB 


- VECTOR CONTAINING 
OF OUTPUT 


MA 


COEFFICIENTS 


FOR 


FORWARD 


PREDICTION 


c 

c 


GA 


- VECTOR CONTAINING 
OF INPUT 


AR 


COEFFICIENTS 


FOR 


BACKWARD 


PREDICTION 


c 

c 


GB 


- VECTOR CONTAINING 
OF INPUT 


MA 


COEFFICIENTS 


FOR 


BACKWARD 


PREDICTION 


c 

c 


ZTA 


- VECTOR CONTAINING 
OF OUTPUT 


AR 


COEFFICIENTS 


FOR 


BACKWARD 


PREDICTION 


c 

c 


ZTB 


- VECTOR CONTAINING 
OF OUTPUT 


MA 


COEFFICIENTS 


FOR 


BACKWARD 


PREDICTION 



C C - ARRAY WHICH STORES CURRENT VALUES OF TXA 

C D - ARRAY WHICH STORES CURRENT VALUES OF TYA 

C E - ARRAY WHICH STORES CURRENT VALUES OF GA 

C F - ARRAY WHICH STORES CURRENT VALUES OF ZTA 

C P - ARRAY WHICH STORES CURRENT VALUES OF TXB 

C Q - ARRAY WHICH STORES CURRENT VALUES OF TYB 

C R - ARRAY WHICH STORES CURRENT VALUES OF GB 

C S - ARRAY WHICH STORES CURRENT VALUES OF ZTB 

C A VECTOR UPDATED AR COEFFICIENTS OF ARMA MODEL 

C B - VECTOR CONTAINING UPDATED MA COEFFICIENTS OF ARMA MODEL. 

C TAU1 - CONSTANT COMPUTED FROM CORRELATION DATA AND PREDICTION 

C ERROR COEFFICIENTS. 

C TAU2 - CONSTANT COMPUTED FROM CORRELATION DATA AND PREDICTION 
C ERROR COEFFICIENTS. 

C TAU3 - CONSTANT COMPUTED FROM CORRELATION DATA AND PREDICTION 
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c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



ERROR COEFFICIENTS. 

TAU4 - CONSTANT COMPUTED FROM CORRELATION DATA AND PREDICTION 
ERROR COEFFICIENTS. 

XMU1 - COEFFICIENT OF AR-TYPE RECURSIVE FORMULA 

YMU1 - COEFFICIENT OF AR-TYPE RECURSIVE FORMULA 

MU 2 - COEFFICIENT OF AR-TYPE RECURSIVE FORMULA 

XMU3 - COEFFICIENT OF AR-TYPE RECURSIVE FORMULA 

YMU3 - COEFFICIENT OF AR-TYPE RECURSIVE FORMULA 

MU4 - COEFFICIENT OF AR-TYPE RECURSIVE FORMULA 

MU5 - COEFFICIENT OF AR-TYPE RECURSIVE FORMULA 

DET - DETERMINANT OF PREDICTION ERROR MATRIX COMPOSED OF 
VX, VY, VXY. 

ERR - ERROR BETWEEN REFERENCE MODEL OUTPUT AND LATTICE 
REALIZATION OUTPUT. 

DIMENSION RX( 0: 2000) ,RY(0: 2000) ,RXY(0: 2000) ,RYX(0: 2000) ,A(0: 21) 
DIMENSION B( 0: 21),TXA(0: 21),TYA(0: 21),TXB(0: 21),TYB(0: 21),GA(0: 21) 
DIMENSION GB( 0: 21) ,ZTA(0: 21) ,ZTB(0: 21) ,C( 0: 21) ,D(0: 21) ,E( 0: 21) 
DIMENSION P( 0: 21) ,Q(0: 21) ,R(0: 21) ,S(0: 21) ,F( 0: 21) 

REAL MU5 ,MU2 ,MU4 



COMPUTE VALUES FOR TAU1 THROUGH TAU4 



T1S=0 

T2S=0 

T3S=0 

T4S=0 

KI=KS-1 

DO 10 1=0, KI 

T1S=T1S-RYX( I+i)*ZTA(KI-I) 
T2S=T2S+RY( I+1)*ZTA(KI-I) 
T3S=T3S+RY( I+1)*GA( I) 
T4S=T4S+RY( 1+1 )*ZTA( I ) 

10 CONTINUE 

T1T=0 
T2T=0 
T3T=0 
T4T=0 

DO 20 J=0 ,KT 

T 1T=T 1T+RX( J+1)*ZTB(KT-J) 
T2T=T2T-RXY( J+l )*ZTB( KT- J ) 
T3T=T3T-RYX(KI -KT+1+J)*GB( J) 
T4T=T4T - RYX( K I -KT+ 1+ J ) * ZTB ( J ) 
20 CONTINUE 

TAU1=T1S+T1T 

TAU2=T2S+T2T 

TAU3=T3S+T3T 

TAU4=T4S+T4T 
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COMPUTE VALUES FOR THE REFLECTION COEFFICIENTS 



C 
C 

XMU1=-TAU1/VZ 

YMU1=-TAU2/VZ 

MU 2 =-VGZ/VZ 

DET =VX*VY-VXY*VXY 

XMU3= - ( VY*TAU 1 - VXY*TAU2 ) / DET 

YMU3= - ( - VXY*TAU 1+VX*TAU2 ) / DET 

MU4 =( VGZ*TAU4-VZ*TAU3 ) / ( VGZ*VGZ -VG*VZ ) 

MU5 =MU4*MU2 
C 

C COMPUTE ARMA COEFFICIENTS 

C 

DO 16 K=0 ,KS 
C(K)=TXA(K) 

D(K)=TYA(K) 

E(K)=GA(K) 

F(K)=ZTA(K) 

16 CONTINUE 

DO 45 J=1 ,KS 
TXA(J)=C(J)+XMU1*F(KS-J) 

TYA ( J ) =D( J ) +YMU 1*F ( KS - J ) 

GA( J ) =E ( J - 1 ) +MU 2*F( J - 1 ) 

ZTA ( J ) =F( J ) +XMU3*C ( KS - J ) +YMU3*D ( KS - J ) +MU4*E ( J - 1 ) +MU5*F( J - 1 ) 
45 CONTINUE 

DO 31 K=0 , KT 
P(K)=TXB(K) 

Q(K)=TYB(K) 

S(K)=ZTB(K) 

R(K)=GB(K) 

31 CONTINUE 

DO 55 J=1,KT 

TXB(J)=P(J)+XMU1*S(KT+1-J) 

TYB( J)=Q( J)+YMU1*S(KT+1- J) 

GB( J) =R( J)+MU2*S( J) 

ZTB(J)=S(J+1)+XMU3*P(KT-J)+YMU3*Q(KT-J)+MU4*R(J)+MU5*S(J) 
55 CONTINUE 

C WRITE (*,176) KS 

C WRITE (9,176) KS 

176 FORMAT ( 12) 

C WRITE (*,175) ( ZTB(K) , K=0,KT) 

C WRITE (9,175) ( ZTB(K) , K=0,KT) 

175 FORMAT (4( 1X,F10. 5) ) 

C 

C UPDATE ERRORS 

C 

650 FORMAT (/' S UPDATE ERROR IS: ' ,F15. 10) 

VX =VX+XMU 1*TAU 1 
VY =VY+YMU1*TAU2 
VXY =VX Y+XMU 1*T AU2 
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VG =VG+MU2*VGZ 

VZ=VZ+( XMU3*TAU 1+YMU3*TAU2 ) +MU4*TAU3+MU5*TAU4 
VGZ=TAU3+MU2*TAU4 
ERR=VY-(VXY**2)/VX 
WRITE (*,650) ERR 
WRITE (9,650) ERR 
WRITE (*,888) VX,VY,VG,VZ 
888 FORMAT (4( 1X.F10. 6) ) 

COMPUTE MODEL COEFFICIENTS 

DO 65 J=1 ,KS 

A(J)=TYA(J)-TXA(J)*VXY/VX 
65 CONTINUE 

DO 70 J=1 ,KT 

B ( J ) =TYB ( J ) -TXB ( J ) *VXY/VX 
70 CONTINUE 

B(0)=-VXY/VX 

RETURN 

END 
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APPENDIX C. SUBROUTINE FOR MAIN PROGRAM 



c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



SUBROUTINE NEWMA(KS , KT , VX , VY , VXY , VG , VZ , VGZ , TXA , TXB , TYA , TYB , GA , GB , Z 
&TA , ZTB , RX , RY , RXY , RYX , A , B ) 

THIS SUBROUTINE COMPUTES MA PARAMETER VALUES FOR AN ARMA MODEL AS 
THE MA ORDER INCREASES BY ONE. 

VARIABLE DEFINITIONS 



KS 

KT 

RX 

RY 

RXY 

RXY 

TXA 

TXB 

TYA 

TYB 

GA 

GB 

ZTA 

ZTB 

C 

D 

E 

F 

P 

Q 

R 

S 

A 

B 

TAU1P 

TAU2P 

TAU3P 



CURRENT AR ORDER 
CURRENT MA ORDER 

AUTOCORRELATION DATA OF INPUT VN 
AUTOCORRELATION DATA OF OUTPUT Y 
CROSS CORRELATION DATA OF INPUT AND OUTPUT 
CROSSCORRELATION DATA OF OUTPUT AND INPUT 
VECTOR CONTAIN! 

OF INPUT. 

VECTOR CONTAINI 
OF INPUT 
VECTOR CONTAINI 
OF OUTPUT 



OF OUTPUT 



OF INPUT 



OF INPUT 



OF OUTPUT 



OF OUTPUT 



AR COEFFICIENTS 


FOR 


FORWARD 


PREDICTION 


MA COEFFICIENTS 


FOR 


FORWARD 


PREDICTION 


AR COEFFICIENTS 


FOR 


FORWARD 


PREDICTION 


MA COEFFICIENTS 


FOR 


FORWARD 


PREDICTION 


AR COEFFICIENTS 


FOR 


BACKWARD 


PREDICTION 


MA COEFFICIENTS 


FOR 


BACKWARD 


PREDICTION 


AR COEFFICIENTS 


FOR 


BACKWARD 


PREDICTION 


MA COEFFICIENTS 


FOR 


BACKYARD 


PREDICTION 


; CURRENT VALUES 


OF 


TXA 




! CURRENT VALUES 


OF 


TYA 




; CURRENT VALUES 


OF 


GA 




: CURRENT VALUES 


OF 


ZTA 




! CURRENT VALUES 


OF 


TXB 




; CURRENT VALUES 


OF 


TYB 




: CURRENT VALUES 


OF 


GB 




; CURRENT VALUES 


OF 


ZTB 




COEFFICIENTS OF 


ARMA MODEL 





VECTOR CONTAINING UPDATED MA COEFFICIENTS OF ARMA MODEL. 
CONSTANT COMPUTED FROM CORRELATION DATA AND PREDICTION 
ERROR COEFFICIENTS. 

CONSTANT COMPUTED FROM CORRELATION DATA AND PREDICTION 
ERROR COEFFICIENTS. 

CONSTANT COMPUTED FROM CORRELATION DATA AND PREDICTION 



6U 



OOU tH C\l OOO 



c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 

c 



ERROR COEFFICIENTS. 

TAU4P - CONSTANT COMPUTED FROM CORRELATION DATA AND PREDICTION 
ERROR COEFFICIENTS. 

XETA1 - COEFFICIENT OF MA-TYPE RECURSIVE FORMULA 

YETA1 - COEFFICIENT OF MA-TYPE RECURSIVE FORMULA 

ETA2 - COEFFICIENT OF MA-TYPE RECURSIVE FORMULA 

XETA3 - COEFFICIENT OF MA-TYPE RECURSIVE FORMULA 

YETA3 - COEFFICIENT OF MA-TYPE RECURSIVE FORMULA 

ETA4 - COEFFICIENT OF MA-TYPE RECURSIVE FORMULA 

ETAS - COEFFICIENT OF MA-TYPE RECURSIVE FORMULA 

DET - DETERMINANT OF PREDICTION ERROR MATRIX COMPOSED OF 
VX,VY,VXY. 

ERR - ERROR BETWEEN REFERENCE MODEL OUTPUT AND LATTICE 

DIMENSION RX( 0: 2000) ,RY(0: 2000) ,RXY(0: 2000) ,RYX(0: 2000) ,A(0: 21) 
DIMENSION B( 0: 21),TXA(0: 21),TXB(0: 21),TYA(0: 21),TYB(0: 21),GA(0: 21) 
DIMENSION GB( 0: 21),ZTA(0: 21),ZTB(0: 21),C(0: 21),D(0: 21) 

DIMENSION E ( 0 : 21),F(0: 21),P(0: 21),Q(0: 21),R(0: 21),S(0: 21) 



COMPUTE VALUES FOR TAU1 PRIME THROUGH TAU4 PRIME 



T1TP=0 

T2TP=0 

T3TP=0 

T4TP=0 

KJ=KT-1 

DO 10 1=0, KJ 

T 1 TP=T 1 TP+RX ( 1+1 )*GB( KJ- I ) 

T2TP=T2TP-RXY( I+1)*GB(KJ-I) 

T3TP=T3TP+RX( I+1)*ZTB( I ) 

T4TP=T4TP+RX( I+1)*GB( I) 

CONTINUE 

T1SP=0 

T2SP=0 

T3SP=0 

T4SP=0 

DO 20 J=0 ,KS 

T1SP=T1SP-RYX(J+1)*GA(KS-J) 

T2SP=T2SP+RY( J+l )*GA( KS - J) 
T3SP=T3SP-RXY(KJ-KS+1+J)*ZTA( J) 
T4SP=T4SP-RXY(KJ-KS+1+J)*GA(J) 

CONTINUE 

TAU1P=T1TP+T1SP 

TAU2P=T2TP+T2SP 

TAU3P=T3TP+T3SP 

TAU4P=T4TP+T4SP 

COMPUTE VALUES FOR REFLECTION COEFFICIENTS 



XETA1=-TAU1P/VG 
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c 

c 

c 



16 

165 



45 



41 

175 



55 

C 

C 

C 



66 



YETA1=-TAU2P/VG 

ETA2 =-VGZ/VG 

DET =VX*VY - VXY*VXY 

XETA 3= - ( V Y -• TAU 1 P - VX Y*TAU 2 P ) / DET 

YETA3=-( -VXY*TAU1P+VX*TAU2P) /DET 

ETA4 =( VGZ*TAU4P-VG*TAU3P ) / ( VGZ*VGZ -VG*VZ ) 

ETA5 =ETA4*ETA2 



COMPUTE ARMA PARAMETERS 

DO 16 K=0 ,KT 
P(K)=TXB(K) 

Q(K)=TYB(K) 

R(K)=GB(K) 

S(K)=ZTB(K) 

CONTINUE 

FORMAT (4( 1X,F10. 5) ) 

DO 45 J=1 ,KT 

TXB(J)=P(J)+XETA1*R(KT-J) 

TYB ( J ) =Q( J ) +YETA 1*R( KT - J ) 

GB( J)=R( J)+XETA3*P(KT-J)+YETA3*Q(KT- J)+ETA4*S( J- 1 )+ETA5*R( J- 1) 
ZTB(J)=S(J-1)+ETA2*R(J-1) 

CONTINUE 
DO 41 K=0 ,KS 
C(K)=TXA(K) 

D(K)=TYA(K) 

E(K)=GA(K) 

F(K)=ZTA(K) 

CONTINUE 

FORMAT (5(1X,F10. 5)) 

DO 55 J=1 ,KS 

TXA(J)=C(J)+XETA1*E(KS+1-J) 

TYA( J)=D( J)+YETA 1*E( KS+1 - J) 

GA( J)=E( J+1)+XETA3*C(KS-J)+YETA3*D(KS-J)+ETA4*F( J)+ETA5*E( J) 
ZTA(J)=F( J)+ETA2*E(J) 

CONTINUE 

UPDATE ERRORS 

VX=VX+XETA1*TAU1P 
VY =VY +YET A 1 *T AU 2 P 
VXY=VXY+XETA 1*TAU 2 P 

VG=VG+(TAU1P*XETA3+TAU2P*YETA3)+ETA4*TAU3P+ETA5*TAU4P 

VZ=VZ+ETA2*VGZ 

VGZ=TAU3P+ETA2*TAU4P 

ERR=VY-(VXY**2)/VX 

WRITE (*,66) ERR 

WRITE (9,66) ERR 

FORMAT (/' T UPDATE ERROR IS: ' ,F15. 10) 
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WRITE (*,889) VX,VY,VG,VZ 
889 FORMAT (4(1X,F10. 6)) 

COMPUTE MODEL COEFFICIENTS 

DO 65 J=1 ,KT 

B(J)=TYB(J)-TXB(J)*VXY/VX 
65 CONTINUE 

DO 70 J=1 ,KS 

A(J)=TYA(J)-TXA(J)*VXY/VX 
70 CONTINUE 

B(0)=-VXY/VX 

RETURN 

END 
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APPENDIX D. ADAPTIVE LATTICE ALGORITHM PROGRAM 



C THIS PROGRAM CALCULATES VALUES OF THE LATTICE COEFFICENTSAND 
C OUTPUT OF AN ARMA DIGITAL LATTICE FILTER USING AN ADAPTIVE 

C LATTICE ALGORITHM. 

C 

C VARIABLE DEFINITIONS 

C 

C X 

C 

C AK 

C R 

C BO 

C M 

C 

C EXF 

C EXB 

C EXBD 

C EYF 

C EYB 

C EYBD 

C ERROR 

C 

C YE 

C 

C MU 

C RHO 

C 

C SIGK 

C 

C SIGR 

C 

C SIGB 

C 
C 

DIMENSION EXF( 10) ,EXB( 10) ,EXBD( 10) ,EYF( 10) ,EYB( 10) ,EYBD( 10) ,R( 10) 
DIMENSION AK(10), X(9900), ERR0R(9900), V(12), YE(6,9900) 

REAL MU 
SIGK=1. 

SIGR=1. 

SIGB=1. 

C 

INITIALIZE ARRAYS 

DO 5 1=1,10 
EXF( I )=0 
EXB( I )=0 



- ARRAY OF INPUT DATA, COMPUTER GENERATED WHITE GAUSSIAN 
NOISE WITH UNIT VARIANCE. 

- ARRAY OF LATTICE COEFFICIENTS. 

- ARRAY OF LATTICE COEFFICIENTS. 

- TERMINAL CONDITION OF LATTICE REALIZATION. 

- NUMBER OF LATTICE STAGES (EQUIVALENT TO ORDER OF ARMA 
MODEL). 

- ARRAY OF FORWARD PREDICTION ERRORS FOR INPUT X. 

- ARRAY OF BACKWARD PREDICTION ERRORS FOR INPUT X. 

- ARRAY OF DELAYED BACKWARD PREDICTION ERRORS FOR INPUT X. 

- ARRAY OF FORWARD PREDICTION ERRORS FOR OUTPUT Y. 

- ARRAY OF BACKWARD PREDICTION ERRORS FOR OUTPUT Y. 

- ARRAY OF DELAYED BACKWARD PREDICTION ERRORS FOR OUTPUT Y. 

- DIFFERENCE BETWEEN REFERENCE MODEL OUTPUT AND LATTICE 
REALIZATION OUTPUT. 

- ARRAYS CONTAINING LATTICE COEFFICIENT VALUES AT EACH 
ITERATION. 

- CONVERGENCE CONSTANT. 

- WEIGHT GIVEN TO CUURENT POWER LEVEL AT EACH STAGE OF THE 
LATTICE STRUCTURE. 

- POWER LEVEL USED TO NORMALIZE CONVERGENCE CONSTANT WHEN 
UPDATING AK LATTICE COEFFICIENTS. 

- POWER LEVEL USED TO NORMALIZE CONVERGENCE CONSTANT WHEN 
UPDATING R LATTICE COEFFICIENTS. 

- POWER LEVEL USED TO NORMALIZE CONVERGENCE CONSTANT WHEN 
UPDATING TERMINAL CONDITION BO. 
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O O O O O O ON 



EXBD( I)=0 
EYF( I)=0 
EYB( I)=0 
EYBD( I )=0 
R(I)=0 
AK(I)=0 

5 CONTINUE 

DO 6 1=1,9000 
X(I)=0 
ERROR( I)=0 
CONTINUE 

ENTER VALUE OF THE CONVERGENCE CONSTANT MU AND VALUE OF RHO. 



M=2 
N=300 

WRITE (6,*) 'ENTER MU' 

RE AD (6,*) MU 
WRITE (6,*) 'ENTER RHO' 

READ (6,*) RHO 

GENERATE WHITE GAUSSIAN NOISE 

ISIZE = N 
IX = 152255 
I SORT = 0 
MUL = 2 

DO 7 K= 1, ISIZE 
CALL SRND( IX, V, 12 ,MUL, ISORT) 

XT=-6. 0 

DO 8 1=1,12 
8 XT=XT+V( I ) 

X(K)=XT 
7 CONTINUE 

C 

C COMPUTE OUTPUT OF REFERENCE MODEL FILTER AND LATTICE STRUCTURE 

C THEN COMPUTE THE ERROR. 

C 

C REFERENCE MODEL 



Y3=0 

Y2=0 

Y1=0 

X3=0 

X2=0 

X1=0 

B0=1. 

DO 100 1=1, N 

C YF=X( I) -0. 8*X1+1. 78*X2+0. 89*Yl-0. 25*Y2 
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YF=0. 5*X( I) -0. 4*X1+. 89*X2+0. 89*Yl-0. 25*Y2 
C YF=X(I)-2. 7*Xl+3. 21*X2-1. 595*X3+1. 95*Y1-1. 62*Y2+0. 54*Y3 
C YF=0. 5*X( I) -0. 95*X1+1. 33*X2-0. 979*X3+1. 69*Yl-0. 962*Y2+0. 2*Y3 
C YF=X(I)+0. 2*Xl-0. 35*X2+1.4*Yl-0. 85*Y2 
C YF=0. 5*X( I) -0. 2*Xl+0. 445*X2+1. 0*Yl-0. 94*Y2 

C Y3=Y2 

Y2=Y1 
Y1=YF 
C X3=X2 

X2=X1 
X1=X(I) 

C 

C LATTICE FILTER 

C 

EXF( 1)=X( I ) 

EXB( 1 )=X( I ) 

DO 10 K=1 ,M 

10 EXF(K+1)=EXF(K)+R(K)*EXBD(K) -AK(K)*EYBD(K) 

EYF(M+1 )=B0*EXF(M+1) 

DO 20 K=1,M 

20 EYF(M+l-K)=EYF(M+2-K)+R(M+l-K)*EXBD(M+l-K)-AK(M+l-K)*EYBD(M+l-K) 

EYB(1)=EYF(1) 

DO 30 K=1,M-1 

EXB ( K+ 1 ) =EXBD ( K ) +R ( K ) *EXF( K ) -R( K ) *E YF ( K ) 

30 EYB(K+1)=EYBD(K)+AK(K)*EYF(K) -AK(K)*EXF(K) 

YL=EYF( 1) 

ERROR( I)=YF-YL 
ERR=YF-YL 

CALL UPDATE ( R , AK , E YBD , EXBD , ERR , MU , M , S IGK , S IGR , BO , RHO) 

CSB=EXF ( M+ 1 ) *EXF ( M+ 1 ) 

SIGB=RHO*SIGB+( l-RHO)*CSB 
B0=B0+(MU/SIGB)*ERR*EXF(M+1) 

DO 40 K=1,M 
EXBD(K)=EXB(K) 

40 EYBD(K)=EYB(K) 

DO 50 J=1,M 
YE( J, I)=AK( J) 

50 YE( J+M, I)=R( J) 

202 FORMAT (2C1X.F10. 6)) 

100 CONTINUE 

C 

C PRINT THE ERROR AND VALUES OF THE LATTICE COEFFICIENTS. 

C 

WRITE (*,200) (ERROR(K) , K=1,N,10) 

WRITE (9,200) (ERROR(K) , K=1,N,10) 

200 FORMAT (5( 1X,F10. 6)) 

WRITE (9,209) 

209 FORMAT( ' ' ) 

WRITE (*,201) (R(K) , K=1,M) 
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205 

201 

C 

C 

C 

C 



C 

C 

C 



20 



10 



C 

C 

C 



10 

C 

C 



WRITE (*,201) ( AK(K) , K=1,M) 

WRITE (9,201) (R(K) , K=1,M) 

WRITE (9,201) (AK(K), K=1,M) 

WRITE (*,205) BO 
WRITE (9,205) BO 
FORMAT (F10.6) 

FORMAT (5( 1X,F10. 6)) 

CALL PLOTTING ROUTINES TO PLOT ERROR AND LATTICE COEFFICIENTS. 

CALL PLOT ( ERROR, N) 

CALL PLOT1 ( YE ,N) 

STOP 

END 

SUBROUTINE WHICH UPDATES LATTICE COEFFICIENTS. 

SUBROUTINE UPDATE ( R , AK , EYBD , EXBD , ERR , MU , M , S IGK , S IGR , BO , RHO ) 
DIMENSION R( 10) ,AK( 10) ,EYBD( 10) ,EXBD( 10) 

REAL MU 
CSK=0. 

CSR=0. 

DO 20 J=1,M 

C SK=C SK+E YBD ( J ) *E YBD ( J ) *B 0** 2 
CSR=CSR+EXBD(J)*EXBD(J)*B0**2 
SIGK=RHO*SIGK+( l-RHO)*CSK 
S IGR=RHO*S IGR+( 1 -RHO) *CSR 
DO 10 J=1,M 

R(J)=R(J)+(MU/SIGR)*ERR*EXBD(J)*BO 

AK( J ) =AK ( J ) - ( MU/S I GK ) *ERR*EYBD( J ) *B 0 

CONTINUE 

RETURN 

END 

PLOTTING ROUTINE TO PLOT ERROR 

SUBROUTINE PLOT(Y,N) 

DIMENSION Y(N) ,X( 9900) 

DO 10 J=1 ,N 
X(J)=J 
CALL TEK618 
CALL PRTPLT( 72,6) 

CALL SHERPAC ADAPTIVE A' ,3) 

CALL RESET(’ALL') 

CALL PAGE ( 8 . 50,6. 0) 

CALL HWROT( 'AUTO' ) 

CALL XINTAX 

CALL AREA2D(5. 0,3. 0) 

CALL HEIGHT(0. 14) 

CALL COMPLX 
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CALL SHDCHR( 90. 0,1,0. 002,1) 

CALL HEADIN( 'LEARNING CURVES? ’ , 100 , 2. 0 , 1) 

CALL XNAME( ' ITERATIONS? ' , 100) 

CALL YNAME( 'ERROR? ' , 100) 

C CALLMESSAGC ADAPTIVE FILTER ?' ,100,3. 0,-0. 8) 

CALL THKFRM(0. 03) 

CALL FRAME 

CALL GRAF (0, 'SCALE' ,N,-3. 00, 'SCALE' ,3. 00) 

C CALL THKCRV(0. 02) 

CALL CURVE(X,Y,N,0) 

CALL ENDPL(O) 

CALL DONEPL 
RETURN 
END 

SUBROUTINE TO PLOT LATTICE COEFFICIENTS 
SUBROUTINE PLOT1 ( YE,N) 

DIMENSION YE(6,9900),X(9900) ,Y(9900) ,YD(9900) ,A( 10) 
TRUE VALUES OF THE PARAMETERS 
A( l)=-0. 240719 
A(2)=0. 125 
A(3)=-0. 195719 
A(4)=0. 8900 
A(5)=0. 89 
DO 10 J=1 ,N 
X(J)=J 
CALL TEK618 
CALL PRTPLT( 72,6) 

CALL SHERPA( 'MENNECKE' ,'A' ,3) 

PRINT SHERPA FILE: SHERPA XXYYZZXX SHGRAPH A 
CALL RESET( 'ALL' ) 

CALL PAGE( 8. 50,6. 0) 

CALL HWROT( 'AUTO' ) 

CALL X INTAX 
CALL AREA2D(5. 0,3. 0) 

CALL HEIGHT(0. 14) 

CALL COMP LX 

CALL SHDCHRC 90. 0,1,0. 002,1) 

CALL HEADINC 'PARAMETERS?' ,100,2. 0,1) 

CALL XNAME( 'ITERATIONS?' ,100) 

CALL YNAME( ' MAGNITUDE? ' , 100) 

CALL MESSAG( ' ADAPTIVE ARMA LATTICE? ' , 100, 3. 0 , -0. 8) 
CALL THKFRM( 0. 03) 

CALL FRAME 

CALL GRAF (0, 'SCALE' ,N,-1. 00, 'SCALE' ,1. 0) 

CALL THKCRV( 0. 02) 

TO PLOT ESTIMATES 
DO 20 K=1 ,4 



6S 



DO 30 J=1,N 
30 Y(J)=YE(K,J) 

CALL CURVE ( X, Y,N, 0) 

20 CONTINUE 

C TO PLOT TRUE PARAMETERS 

DO 40 K=1 ,4 
DO 50 J=1,N 
50 Y(J)=A(K) 

CALL DASH 

CALL CURVE (X,Y,N,0) 

40 CONTINUE 

CALL ENDPL(O) 

CALL DONE PL 

RETURN 

END 
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APPENDIX E. DERIVATION OF DIFFERENCE EQUATION FOR TWO 

STAGE ARMA LATTICE FILTER 

The lattice filter is described by the expressions for the forward and backward pre- 
diction errors, equations (3.6) and (3.7) respectively, 



ey (A) = x(k) + iv 2 x(k - 1) - vv, 1 y{k - 1) 

% ( k ) = % (*) + w 2 % ( k ~ 0 - lv i ( k ~ 1) 

% i k ) = x(k - I) + iv -2 x(k) - vvj y(k) (E - 1) 

4, W =y( k - J ) - M k ) + »3 }'( k ) 

e A ( A ') = € fi W + w l % ( k ~ 0 - u 3 % ( k ~ 0 

and the expression for the filter output. 

y(k) = ey (A) + iv 4 x(k - 1) - vvj y(k - 1) {E - 2) 



Substituting for (A') in equation (E-2) yields. 



y(A) = e* (A) + vv 4 (A - 1) - iv 3 2 e } y (A - I) + vv 4 x{k - 1) - vvj y{k - 1) (£ - 3) 

Substituting for eg (A — 1) and eg x (A — 1) in equation (E-3) 

r(A) = e'} 2 (A) + vv 2 [ .v(A - 2) + iv 1 x(k - 1) - y(k - I) ] 

~ lv 3 [ y( k - 2) - wj x(k - 1) -f wj y(k - I) ] (E - 4) 

+ x(k - 1) - vv] v(A - 1) 



Substituting for e} 2 (A) in (E-4) we obtain 

><A) = Cy (A) + w\ c'y (A - 1) - w 2 e b] (A - 1) + w 4 x(k - 2) + vv 4 2 ns x(k - 1) 

—iv 2 vv] y(k — 1) — IV 2 y{k — 2) + vv 2 w] x(k — 1) — vv] vv] y(k — 1) (£— 5) 

+vvj x(k — I) — vv] y(k — 1) 

From (E-l). substitute for el (A — I) and eg (A — 1) in (E-5) to obtain 
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/ 



(£- 6 ) 



>•(*) = % ( k ) + w 2 [ x ( k ~-) + w 2 x ( k ~ 1 ) ~ w l A k - 1 ) ] 

- w, 2 [ y{k - 2 ) - wj x(k - 1) + wj y(k - 1) ] 
4- vi'j x(k — 2 ) 4- w 2 ivs x(k — 1 ) — w? iv! y(k — 1 ) 
-wj y(k - 2) + vv| ivj x(k - 1) - w 3 2 w' 2 y(k - 1) 

+ wl x(k - 1) - U3 1 y(k - 1) 



Now, from (E-l), substituting for e) {k) in (E-6) yields 

y(A) = x(k) + vv - 2 x(k — 1) — w' y(k — 1) + w 2 x(k — 2) + ivj tvj x(k — 1) 

— wl iv 4 ' y(k — 1) —wly(k — 2) + wf w' x(k — 1) — w 2 w\ y(k — 1) 

+ iv? x(k — 2 ) 4- wj u-2 x(k — 1) —wj vvj y{k — 1) — iv 2 y(k — 2 ) 

4- U-? wl x(k — 1) — vv 3 u-j y(k — 1)4- vvj x(k — 1) — u-] y(k — 1 ) 



Grouping terms we get, 



y(k) = x(k) 4- (ivj 4- u -2 wl + vv’ vv 2 + wj wj + wl vv 2 4 - vv]) x(k — 1) 

+ ( u 2 + w b x ( k - 2) 

/ 1 , 12, 12. 12. 12, U n i\ 

— (vv, + U 4 u -2 4- vv 3 IV, 4- iv 4 vv 4 4- vv 3 vv 3 4- iv 3 ) v v(A: — 1) 

- (vv 2 + w 2 3 )y(k - 2) 



(E- S) 



From the gradient estimator and coefficient update equations we know that the follow- 
ing relationships among lattice coefficients are true 



1 1 

H'2 = Uj 



2 

v 4 



(£- 9 ) 



Using these equalities in equation (E-8). we obtain the final expression for the difference 
equation. 



y{k) = x (k) + 2 (h -2 + wl wl + wj wj) x(k - 1) + 2 vv 2 2 x(k - 2) 
— 2(vv, -)- iv *2 w '2 4- vv'-| vvj* j y(k — 1) — 2 vv, y\k — 2) 



(£— 10) 
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